      SUBROUTINE DVDQ (NEQ,T,Y,F,KD,EP,IFLAG,H,HMINA,HMAXA,DELT,        DVDQ0001
     1  TFINAL,MXSTEP,KSTEP,KEMAX,EMAX,KQ,YN,DT)                        DVDQ0002
C                                                                       DVDQ0003
C     DOUBLE PRECISION VARIABLE ORDER INTEGRATION SUBROUTINE            DVDQ0004
C     FOR THE SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS.              DVDQ0005
C                                                                       DVDQ0006
C     ANALYSIS AND CODING BY FRED T. KROGH,  AT THE JET PROPULSION      DVDQ0007
C     LABORATORY, PASADENA, CALIF.  APRIL 1, 1969.                      DVDQ0008
C     THIS SUBROUTINE DESIGNED FOR THE UNIVAC 1108.                     DVDQ0009
C                                                                       DVDQ0010
C     AT THE END OF THIS LISTING INSTRUCTIONS ARE GIVEN FOR REMOVING    DVDQ0011
C     SOME FEATURES AND FOR ADDING OTHERS. THE GSTOP FEATURE IS         DVDQ0012
C     EXPLAINED NEAR THE END OF THE LISTING.                            DVDQ0013
C                                                                       DVDQ0014
C     VARIABLES IN THE CALLING SEQUENCE HAVE THE FOLLOWING TYPES.       DVDQ0015
      INTEGER NEQ,KD(1),IFLAG,MXSTEP,KSTEP,KEMAX,KQ(6)                  DVDQ0016
      REAL    EP(6),HMINA,HMAXA,EMAX                                    DVDQ0017
      DOUBLE PRECISION T,Y(6),F(6),H,DELT,TFINAL,YN(6),DT(20,6)         DVDQ0018
C                                                                       DVDQ0019
C     PARAMETERS WHICH MUST BE ASSIGNED VALUES BEFORE CALLING           DVDQ0020
C     DVDQ ARE    NEQ, T, Y, KD, H, HMINA, HMAXA, DELT,                 DVDQ0021
C                 TFINAL, AND MXSTEP.                                   DVDQ0022
C     DVDQ IS USED ONLY ON THE INITIAL ENTRY. ALL OTHER                 DVDQ0023
C     ENTRIES ARE MADE BY CALLING DVDQ1. IN ADDITION TO                 DVDQ0024
C     THE PARAMETERS MENTIONED ABOVE THE USER MUST ASSIGN               DVDQ0025
C     VALUES TO F (ONCE PER STEP INITIALLY, AND TWICE PER STEP          DVDQ0026
C     AFTER GETTING STARTED) AND EP (EITHER INITIALLY, OR DURING        DVDQ0027
C     THE INTEGRATION IF A RELATIVE ERROR TEST IS USED).                DVDQ0028
C     THE FOLLOWING PARAMETERS GIVE ADDITIONAL INFORMATION ABOUT THE    DVDQ0029
C     INTEGRATION AND ARE USED FOR STORAGE. THEY SHOULD NOT BE          DVDQ0030
C     CHANGED BY THE USER.  IFLAG,KSTEP,KEMAX,EMAX,KQ,YN, AND DT.       DVDQ0031
C                                                                       DVDQ0032
C                                                                       DVDQ0033
C                                                                       DVDQ0034
C     AN EXAMPLE OF HOW ONE MIGHT SET UP THE CALLS TO DVDQ IS GIVEN     DVDQ0035
C     BELOW.                                                            DVDQ0036
C     CALL DVDQ(NEQ,...,DT)                                             DVDQ0037
C     GO TO YYY                                                         DVDQ0038
C XXX CALL DVDQ1                                                        DVDQ0039
C YYY GO TO (N1,N2,...,N8), IFLAG                                       DVDQ0040
C     AFTER THE INITIAL CALL, RETURN TO DVDQ WITH                       DVDQ0041
C     GO TO XXX                                                         DVDQ0042
C     IF NO ERROR HAS BEEN MADE, IFLAG WILL EQUAL 1 AFTER THE INITIAL   DVDQ0043
C     CALL INDICATING THE DERIVATIVES ARE TO BE COMPUTED.               DVDQ0044
C                                                                       DVDQ0045
C                                                                       DVDQ0046
C     THE USAGE OF THE VARIABLES IS GIVEN BELOW.                        DVDQ0047
C                                                                       DVDQ0048
C     NEQ=NUMBER OF EQUATIONS (INPUT)                                   DVDQ0049
C                                                                       DVDQ0050
C     T=INDEPENDENT VARIABLE (INITIAL VALUE SUPPLIED BY THE USER)       DVDQ0051
C                                                                       DVDQ0052
C     Y=CURRENT VALUE OF DEPENDENT VARIABLE. THE INITIAL                DVDQ0053
C       VALUE OF Y MUST BE SPECIFIED BY THE USER BEFORE                 DVDQ0054
C       THE FIRST ENTRY. THE DIMENSION OF Y MUST BE                     DVDQ0055
C       AT LEAST AS GREAT AS THE SUM OF THE ORDERS OF                   DVDQ0056
C       THE DIFFERENTIAL EQUATIONS WHICH ARE BEING                      DVDQ0057
C       INTEGRATED. IF WE LET KD(I) DENOTE THE ORDER                    DVDQ0058
C       OF THE I-TH DIFFERENTIAL EQUATION, THEN Y(J)                    DVDQ0059
C       IS THE K-TH DERIVATIVE OF THE L-TH COMPONENT,                   DVDQ0060
C       WHERE L IS THE SMALLEST INTEGER FOR WHICH                       DVDQ0061
C       KD(1)+KD(2)+...+KD(L).GE.J  AND K=KD(L)+J-1-(KD(1)              DVDQ0062
C       +KD(2)+...+KD(L)), J=1,2,...,(KD(1)+KD(2)+...+KD(NEQ)).         DVDQ0063
C       (FOR EXAMPLE, FOR THE SYSTEM  F(1)=UPP, F(2)=VPP, WHERE P       DVDQ0064
C       DENOTES A PRIME, Y(1)=U, Y(2)=UP, Y(3)=V, Y(4)=VP.)             DVDQ0065
C                                                                       DVDQ0066
C     F(I)=KD(I)-TH DERIVATIVE OF THE I-TH COMPONENT WITH RESPECT       DVDQ0067
C       TO T, I=1,2,...,NEQ. THE USER MUST PROVIDE                      DVDQ0068
C       THE CODE WHICH COMPUTES F GIVEN Y AND T.                        DVDQ0069
C                                                                       DVDQ0070
C     KD GIVES THE ORDER OF THE DIFFERENTIAL EQUATIONS IN THE           DVDQ0071
C        SYSTEM. KD MUST BE LESS THAN OR EQUAL TO 4.                    DVDQ0072
C        (FOR DIFFERENTIAL EQUATIONS WITH DIFFERENT ORDERS SET          DVDQ0073
C        KD.LT.0. IF THIS IS DONE IT IS ASSUMED THAT KD IS A VECTOR     DVDQ0074
C        AND THAT ABS(KD(I)) GIVES THE ORDER OF THE I-TH EQUATION.)     DVDQ0075
C                                                                       DVDQ0076
C     EP IS A PARAMETER USED TO CONTROL THE LOCAL ERROR.                DVDQ0077
C     IF EP IS POSITIVE THE LOCAL ERROR IS KEPT LESS                    DVDQ0078
C     THAN EP IN ALL COMPONENTS OF THE DIFF. EQ.                        DVDQ0079
C     (THE ESTIMATED LOCAL ERROR IS KEPT LESS THAN EP IN                DVDQ0080
C     THE (KD(I)-1)-ST DERIVATIVE OF THE I-TH COMPONENT. THUS           DVDQ0081
C     FOR EQUATIONS WITH ORDER GREATER THAN ONE, THE ERROR              DVDQ0082
C     IN A DERIVATIVE IS ESTIMATED. IN THIS CASE THE VALUE OF           DVDQ0083
C     EP REQUIRED TO OBTAIN A GIVEN ACCURACY IN THE DEPENDENT           DVDQ0084
C     VARIABLE DEPENDS ON THE SCALING.)                                 DVDQ0085
C     IF EP.LT.0, THEN IT IS ASSUMED THAT EP                            DVDQ0086
C     IS A VECTOR. LET K BE THE SMALLEST VALUE                          DVDQ0087
C     OF I FOR WHICH EP(I).GE.0.  FOR I.LT.K                            DVDQ0088
C     THE LOCAL ERROR CONTROL IS BASED ON                               DVDQ0089
C     ABS(EP(I)), AND FOR I.GE.K IT IS BASED ON                         DVDQ0090
C     EP(K). IF ONE WANTS A RELATIVE ERROR TEST--                       DVDQ0091
C     FOR EXAMPLE, THE LOCAL ERROR IS TO BE KEPT                        DVDQ0092
C     LESS THAN C*P WHERE C IS A CONSTANT                               DVDQ0093
C     AND P IS A POSITIVE FUNCTION OF T AND Y,                          DVDQ0094
C     THEN ONE SHOULD SET EP=C*P WHEN IFLAG=1.                          DVDQ0095
C     IF EP=0 AND HMAXA.NE.0, IFLAG IS SET EQUAL 8. IF EP=0 AND HMAXA=0,DVDQ0096
C     NO ERROR TESTS ARE PERFORMED AND THE ORDER(S) AND STEPSIZE ARE    DVDQ0097
C     NOT CHANGED. THIS OPTION SHOULD NOT BE USED IF KQ(I)=1 FOR ANY I. DVDQ0098
C                                                                       DVDQ0099
C     IFLAG IS USED FOR COMMUNICATION BETWEEN THE INTEGRATOR            DVDQ0100
C           AND THE PROGRAM WHICH CALLS IT.  THE VALUE                  DVDQ0101
C           OF IFLAG SHOULD NOT BE CHANGED BY THE USER.                 DVDQ0102
C        THE FOLLOWING VALUES OF IFLAG HAVE THE FOLLOWING MEANINGS.     DVDQ0103
C     =1  THE VALUE OF Y FOR THE CURRENT STEP HAS BEEN                  DVDQ0104
C         PREDICTED. THE USER SHOULD COMPUTE F AND CALL DVDQ1.          DVDQ0105
C         IF A RELATIVE ERROR TEST IS USED THE NEW VALUE                DVDQ0106
C         OF EP SHOULD ALSO BE COMPUTED HERE.                           DVDQ0107
C     =2  THE VALUE OF Y FOR THE CURRENT STEP HAS BEEN                  DVDQ0108
C         CORRECTED. THE USER SHOULD COMPUTE F AND CALL DVDQ1.          DVDQ0109
C     =3  AN OUTPUT POINT HAS BEEN REACHED (SEE DESCRIPTION             DVDQ0110
C         OF DELT), PRINT RESULTS AND CALL DVDQ1.                       DVDQ0111
C     =4  T=TFINAL    IF DVDQ1 IS CALLED WITH T=TFINAL AND              DVDQ0112
C         IFLAG=4, IFLAG IS SET EQUAL TO 8. IF THE VALUE OF             DVDQ0113
C         TFINAL IS CHANGED THE INTEGRATION WILL CONTINUE.              DVDQ0114
C     =5  KSTEP=KSOUT   (SEE THE DESCRIPTION OF MXSTEP).                DVDQ0115
C     =6  EMAX.GT..1 AND IT APPEARS TO THE SUBROUTINE THAT              DVDQ0116
C         REDUCING H WILL NOT HELP BECAUSE OF ROUND-OFF ERROR.          DVDQ0117
C         IF THIS OCCURS A LARGER VALUE OF EP (OR OF ABS(EP(KEMAX)) IF  DVDQ0118
C         EP IS A VECTOR) SHOULD PROBABLY BE USED. IF EP IS NOT         DVDQ0119
C         INCREASED, TOO SMALL A STEPSIZE IS LIABLE TO BE USED. (WE HAVEDVDQ0120
C         FOUND THAT REPLACING EP WITH 32.*EMAX*EP WORKS QUITE WELL.)   DVDQ0121
C         INCREASING EP IN THIS WAY WILL NOT DEGRADE THE ACCURACY,      DVDQ0122
C         HOWEVER IF THE NATURE OF THE PROBLEM CHANGES IT MAY PAY TO    DVDQ0123
C         USE A SMALLER VALUE OF EP LATER IN THE INTEGRATION.           DVDQ0124
C     =7  ABS(H).LT.HMINA.  TO CONTINUE WITH THE CURRENT                DVDQ0125
C         VALUE OF H, SET HMINA.LE.ABS(H) AND CALL DVDQ1.               DVDQ0126
C         IF THE INTEGRATOR HAS JUST HALVED H ONE MAY CONTINUE          DVDQ0127
C         WITH TWICE THE STEPSIZE BY SIMPLY CALLING DVDQ1. (SUCH        DVDQ0128
C         AN ACTION IS RISKY WITHOUT A CAREFUL ANALYSIS OF THE          DVDQ0129
C         SITUATION.) IF THE STEPSIZE HAS NOT JUST BEEN HALVED          DVDQ0130
C         (ABS(H).LT.HMINA MAY BE DUE TO THE USER INCREASING THE        DVDQ0131
C         VALUE OF HMINA OR TO HAVING TOO SMALL AN H AT THE END         DVDQ0132
C         OF THE STARTING PHASE.) THE INTEGRATION WILL CONTINUE         DVDQ0133
C         WITH THE CURRENT VALUE OF H AND A RETURN TO THE USER WITH     DVDQ0134
C         IFLAG=7 WILL BE MADE ON EVERY STEP UNTIL ABS(H).GE.HMINA.     DVDQ0135
C     =8  ILLEGAL PARAMETER IN THE CALLING SEQUENCE. IF DVDQ1           DVDQ0136
C         IS CALLED WITH IFLAG=8 THE PROGRAM IS STOPPED.                DVDQ0137
C                                                                       DVDQ0138
C     H=CURRENT VALUE OF THE STEPSIZE. IN SELECTING THE INITIAL         DVDQ0139
C       VALUE FOR H, THE USER SHOULD REMEMBER THE FOLLOWING--           DVDQ0140
C      1. THE INTEGRATOR IS CAPABLE OF CHANGING H QUITE QUICKLY AND     DVDQ0141
C         THUS THE INITIAL CHOICE IS NOT CRITICAL.                      DVDQ0142
C      2. IF IT DOES NOT LEAD TO PROBLEMS IN COMPUTING THE DERIVATIVES  DVDQ0143
C         (E.G. BECAUSE OF OVERFLOW OR TRYING TO EXTRACT THE SQUARE     DVDQ0144
C         ROOT OF A NEGATIVE NUMBER), IT IS BETTER TO CHOOSE H MUCH     DVDQ0145
C         TOO LARGE THAN MUCH TOO SMALL.                                DVDQ0146
C      3. IF H*DELT.LE.0 INITIALLY, AN IMMEDIATE RETURN IS MADE         DVDQ0147
C         WITH IFLAG=8. THE SIGN OF H IS WHAT DETERMINES THE            DVDQ0148
C         DIRECTION OF INTEGRATION.                                     DVDQ0149
C      4. IF DELT=H*(2**K) K A NONNEGATIVE INTEGER THEN OUTPUT          DVDQ0150
C         VALUES WILL BE OBTAINED WITHOUT DOING AN INTERPOLATION.       DVDQ0151
C                                                                       DVDQ0152
C     HMINA   AFTER GETTING STARTED, AND WHENEVER H                     DVDQ0153
C             IS HALVED, ABS(H) IS COMPARED WITH HMINA.                 DVDQ0154
C             IF ABS(H).LT.HMINA CONTROL IS RETURNED TO                 DVDQ0155
C             THE USER WITH IFLAG=7.                                    DVDQ0156
C                                                                       DVDQ0157
C     HMAXA   THE STEPSIZE IS NOT DOUBLED IF                            DVDQ0158
C             DOING SO WOULD MAKE ABS(H).GT.HMAXA                       DVDQ0159
C                                                                       DVDQ0160
C     DELT  ENABLES THE USER TO SPECIFY THE POINTS WHERE                DVDQ0161
C       OUTPUT IS DESIRED. LET TOUT=DELT + THE VALUE OF T THE LAST      DVDQ0162
C       TIME CONTROL WAS RETURNED TO THE USER WITH IFLAG=3. (INITIALLY  DVDQ0163
C       TOUT=THE INITIAL VALUE OF T.) CONTROL IS RETURNED TO THE        DVDQ0164
C       USER WITH IFLAG=3 WHENEVER T=TOUT. IF TOUT DOES NOT FALL        DVDQ0165
C       ON AN INTEGRATION STEP, OUTPUT VALUES ARE OBTAINED BY           DVDQ0166
C       INTERPOLATION ON THE FIRST STEP THAT (T-TOUT)*H.GT.0.           DVDQ0167
C       INTERPOLATED VALUES FOR BOTH Y AND F ARE COMPUTED.              DVDQ0168
C       (NOTE THAT A RETURN WITH IFLAG=3 IS ALWAYS MADE                 DVDQ0169
C       BEFORE TAKING THE FIRST STEP.)                                  DVDQ0170
C                                                                       DVDQ0171
C     TFINAL   CONTROL IS RETURNED TO THE USER WITH IFLAG=4 WHEN        DVDQ0172
C       T REACHES TFINAL. IF TFINAL DOES NOT FALL ON AN INTEGRATION     DVDQ0173
C       STEP VALUES AT TFINAL ARE OBTAINED BY EXTRAPOLATION.            DVDQ0174
C                                                                       DVDQ0175
C     MXSTEP   ON THE INITIAL ENTRY, AND ON ENTRIES                     DVDQ0176
C       WITH  2.LT.IFLAG.LT.6  KSOUT IS SET EQUAL TO                    DVDQ0177
C       KSTEP+MXSTEP. AT THE END OF EACH STEP KSTEP IS INCREMENTED      DVDQ0178
C       AND COMPARED WITH KSOUT. IF KSTEP.GE.KSOUT CONTROL IS           DVDQ0179
C       RETURNED TO THE USER WITH IFLAG=5. (THUS IF DELT IS             DVDQ0180
C       SUFFICIENTLY LARGE, CONTROL WILL BE RETURNED TO THE USER        DVDQ0181
C       WITH IFLAG=5 EVERY MXSTEP STEPS.)                               DVDQ0182
C                                                                       DVDQ0183
C     KSTEP=NUMBER OF INTEGRATION STEPS TAKEN (COMPUTED                 DVDQ0184
C       BY THE INTEGRATOR.)                                             DVDQ0185
C                                                                       DVDQ0186
C     KEMAX=INDEX OF COMPONENT RESPONSIBLE FOR THE                      DVDQ0187
C           VALUE OF EMAX (SEE BELOW).                                  DVDQ0188
C                                                                       DVDQ0189
C     EMAX=LARGEST VALUE IN ANY COMPONENT OF (ESTIMATED ERROR)/EP       DVDQ0190
C       ORDINARILY THE STEPSIZE IS HALVED IF EMAX.GT..1. WITH A         DVDQ0191
C       RECENT HISTORY OF LOCAL ROUND-OFF PROBLEMS VALUES OF EMAX AS    DVDQ0192
C       LARGE AS 1 ARE PERMITTED. THE STEPSIZE IS NOT HALVED ON ANY     DVDQ0193
C       STEP THAT ROUND OFF ERROR APPEARS TO BE LIMITING THE PRECISION. DVDQ0194
C                                                                       DVDQ0195
C     KQ(I)=HIGHEST ORDER DIFFERENCE USED IN INTEGRATING                DVDQ0196
C       THE I-TH EQUATION. (COMPUTED BY THE INTEGRATOR)                 DVDQ0197
C                                                                       DVDQ0198
C     YN=A VECTOR WITH THE DIMENSION OF Y USED TO STORE                 DVDQ0199
C       THE VALUE OF Y AT THE END OF EACH INTEGRATION STEP.             DVDQ0200
C                                                                       DVDQ0201
C     DT=AN ARRAY WITH DIMENSION DT(20,NEQ) USED TO                     DVDQ0202
C       STORE THE DIFFERENCE TABLE.                                     DVDQ0203
C                                                                       DVDQ0204
C                                                                       DVDQ0205
      DOUBLE PRECISION TOUT,TL,TPD,TPD1,TPD2,HH,FAC                     DVDQ0206
      DOUBLE PRECISION DD,D,GAM,GAS,PT,TP                               DVDQ0207
      DIMENSION DD(22),D(21),GAM(20,4),GAS(20),PT(21),FAC(3)            DVDQ0208
      EQUIVALENCE(DD(2),D(1))                                           DVDQ0209
      DIMENSION ETA(18,18)                                              DVDQ0210
C                                                                       DVDQ0211
C*** DECLARATION OF VARIABLES USED IN ENTRY DVDQG(NG,NSTOP,G,GT) ***
C
      INTEGER NG,NSTOP                                                  DVDQ1193
      DOUBLE PRECISION G(1),GT(1)                                       DVDQ1194
      DOUBLE PRECISION RG                                               DVDQ1264
      DOUBLE PRECISION GI                                               DVDQ1265
      DIMENSION GI(2),RG(3)                                             DVDQ1266
C
C
      DATA KMAXO/4/                                                     DVDQ0212
C     KMAXO IS THE MAXIMUM ORDER DIFFERENTIAL EQUATION                  DVDQ0213
C     THIS SUBROUTINE WILL INTEGRATE.                                   DVDQ0214
C                                                                       DVDQ0215
      DATA FAC/1.D0,.5D0,.1666666666666666667D0/                        DVDQ0216
C     FAC(J)=1/(FACTORIAL J), J=1,2,...,MAX(2,KMAXO-1)                  DVDQ0217
C                                                                       DVDQ0218
      DATA KQMAX/19/                                                    DVDQ0219
C     KQMAX GIVES THE MAXIMUM ORDER.                                    DVDQ0220
C     THERE IS LITTLE POINT IN HAVING KQMAX MUCH BIGGER THAN THE NUMBER DVDQ0221
C     OF DECIMAL DIGITS IN THE MANTISSA.                                DVDQ0222
C     IF KQMAX IS SET LESS THAN 6, DT, D, AND PT SHOULD BE DIMENSIONED  DVDQ0223
C     AS IF KQMAX=6.                                                    DVDQ0224
C                                                                       DVDQ0225
      DATA RND,KBIT2/6.9E-18,122/                                       DVDQ0226
C     RND IS APPROXIMATELY 2**(3-B) WHERE B IS                          DVDQ0227
C     THE NUMBER OF BITS IN THE MANTISSA.                               DVDQ0228
C     KBIT2=2*B+2 WHERE B IS THE NUMBER OF BITS IN THE MANTISSA.        DVDQ0229
C     IF THE DERIVATIVES ARE NOT COMPUTED TO THE ACCURACY EXPECTED      DVDQ0230
C     FROM THE WORD LENGTH OF THE COMPUTER (FOR EXAMPLE BECAUSE OF      DVDQ0231
C     CANCELLATION PROBLEMS OR TABULAR DATA), THEN THESE CONSTANTS      DVDQ0232
C     CAN BE CHANGED TO REFLECT THE NUMBER OF BITS WHICH ARE            DVDQ0233
C     SIGNIFICANT IN THE COMPUTED DERIVATIVES. (THIS IS NOT NECESSARY,  DVDQ0234
C     BUT IS WISE IF THE ACCURACY REQUESTED IS DIFFICULT TO OBTAIN      DVDQ0235
C     BECAUSE THE DERIVATIVES HAVE SO FEW SIGNIFICANT DIGITS.)          DVDQ0236
C                                                                       DVDQ0237
      DATA P1,P01,P25,P3E1/.1,.01,.25,3./                               DVDQ0238
C     THE ABOVE DATA STATEMENT CONTAINS VARIOUS CONSTANTS               DVDQ0239
C     USED IN THE SUBROUTINE.                                           DVDQ0240
C                                                                       DVDQ0241
      DATA PT/1.D0,2.D0,4.D0,8.D0,16.D0,32.D0,64.D0,128.D0,256.D0,      DVDQ0242
     1  512.D0,1024.D0,2048.D0,4096.D0,8192.D0,16384.D0,32768.D0,       DVDQ0243
     2  65536.D0,131072.D0,262144.D0,524288.D0,1048576.D0/              DVDQ0244
C     PT(J)=2**(J-1), J=1,2,...,KQMAX+2                                 DVDQ0245
      DATA PTS1,PTS2,PTS3,PTS4,PTS5,P5/1.,2.,4.,8.,16.,.5/              DVDQ0246
C                                                                       DVDQ0247
      DATA GAS/1.D0,-.5D0,-8.333333333333333333D-02,                    DVDQ0248
     1  -4.16666666666666667D-02,-2.63888888888888889D-02,              DVDQ0249
     2  -1.875D-02,              -1.42691798941798942D-02,              DVDQ0250
     3  -1.13673941798941799D-02,-9.35653659611992945D-03,              DVDQ0251
     4  -7.89255401234567901D-03,-6.78584998463470686D-03,              DVDQ0252
     5  -5.92405641233766234D-03,-5.23669325795028507D-03,              DVDQ0253
     6  -4.67749840704226452D-03,-4.21495223900547286D-03,              DVDQ0254
     7  -3.82689955321188442D-03,-3.49734984534991765D-03,              DVDQ0255
     8  -3.21449643132356745D-03,-2.96944771545820961D-03,              DVDQ0256
     9  -2.75539029943671585D-03/                                       DVDQ0257
C     GAS(I) GIVES THE I-TH ADAMS-MOULTON CORRECTOR                     DVDQ0258
C     COEFFICIENT, I=1,2,...,KQMAX+1.                                   DVDQ0259
C                                                                       DVDQ0260
      DATA (GAM(I,1), I=1,20)/                                          DVDQ0261
     1  1.D0,.5D0,.4166666666666666667D0,.375D0,                        DVDQ0262
     2  .348611111111111111D0,.329861111111111111D0,                    DVDQ0263
     3  .315591931216931217D0,.304224537037037037D0,                    DVDQ0264
     4  .294868000440917108D0,.286975446428571429D0,                    DVDQ0265
     5  .280189596443936722D0,.274265540031599059D0,                    DVDQ0266
     6  .269028846773648774D0,.264351348366606510D0,                    DVDQ0267
     7  .260136396127601037D0,.256309496574389153D0,                    DVDQ0268
     8  .252812146729039235D0,.249597650297715667D0,                    DVDQ0269
     9  .246628202582257458D0,.243872812282820742D0/                    DVDQ0270
      DATA (GAM(I,2), I=1,20)/                                          DVDQ0271
     1  .5D0,.1666666666666666667D0,.125D0,.1055555555555555556D0,      DVDQ0272
     2  9.375D-2,              8.56150793650793651D-2,                  DVDQ0273
     3  7.95717592592592593D-2,7.48522927689594356D-2,                  DVDQ0274
     4  7.10329861111111111D-2,6.78584998463470686D-2,                  DVDQ0275
     5  6.51646205357142857D-2,6.28403190954034208D-2,                  DVDQ0276
     6  6.08074792915494387D-2,5.90093313460766200D-2,                  DVDQ0277
     7  5.74034932981782663D-2,5.59575975255986825D-2,                  DVDQ0278
     8  5.46464393325006467D-2,5.34500588782477730D-2,                  DVDQ0279
     9  5.23524156892976012D-2,5.13404509001447570D-2/                  DVDQ0280
      DATA (GAM(I,3), I=1,20)/                                          DVDQ0281
     1  .1666666666666666667D0,4.16666666666666667D-2,                  DVDQ0282
     2  2.91666666666666667D-2,2.36111111111111111D-2,                  DVDQ0283
     3  2.03373015873015873D-2,1.81299603174603175D-2,                  DVDQ0284
     4  1.65181327160493827D-2,1.52772266313932981D-2,                  DVDQ0285
     5  1.42851881914381914D-2,1.34693965531639143D-2,                  DVDQ0286
     6  1.27836579217097570D-2,1.21970388231238926D-2,                  DVDQ0287
     7  1.16879616455733216D-2,1.12408663352884755D-2,                  DVDQ0288
     8  1.08442182943468791D-2,1.04892655447842863D-2,                  DVDQ0289
     9  1.01692338611494262D-2,9.87878870055154653D-3,                  DVDQ0290
     $  9.61366549695202001D-3,9.37040949739634796D-3/                  DVDQ0291
      DATA (GAM(I,4), I=1,20)/                                          DVDQ0292
     1  4.16666666666666667D-2,8.33333333333333333D-3,                  DVDQ0293
     2  5.55555555555555556D-3,4.36507936507936508D-3,                  DVDQ0294
     3  3.67890211640211640D-3,3.22365520282186949D-3,                  DVDQ0295
     4  2.89544753086419753D-3,2.64543583988028432D-3,                  DVDQ0296
     5  2.44737491482283149D-3,2.28579543818052416D-3,                  DVDQ0297
     6  2.15093669481483635D-3,2.03630871020228384D-3,                  DVDQ0298
     7  1.93741301123433302D-3,1.85102419106078320D-3,                  DVDQ0299
     8  1.77476374781296400D-3,1.70683564605258723D-3,                  DVDQ0300
     9  1.64585591005465158D-3,1.59073922159715914D-3,                  DVDQ0301
     $  1.54062133051925634D-3,1.49480466961076253D-3/                  DVDQ0302
C     GAM(I,J) GIVES THE I-TH ADAMS-FALKNER PREDICTOR                   DVDQ0303
C     COEFFICIENT FOR INTEGRATING J-TH ORDER DIFFERENTIAL               DVDQ0304
C     EQUATIONS, I=1,2,...,KQMAX+1,  J=1,2,...,KMAXO.                   DVDQ0305
C                                                                       DVDQ0306
      DATA (ETA(I, 1), I=1,18)/  3.33333330E-01, 2.50000000E-01,        DVDQ0307
     1  1.13636360E-01, 6.73076930E-02, 4.60526330E-02, 3.43749980E-02, DVDQ0308
     2  2.71381590E-02, 2.22547310E-02, 1.87484580E-02, 1.61123220E-02, DVDQ0309
     3  1.40603000E-02, 1.24197060E-02, 1.10802170E-02, 9.96793590E-03, DVDQ0310
     4  9.03137260E-03, 8.23348160E-03, 7.54686710E-03, 6.95082750E-03/ DVDQ0311
      DATA (ETA(I, 2), I=1,18)/  2.00000000E-01, 4.00000000E-01,        DVDQ0312
     1  3.40909090E-01, 2.01923080E-01, 1.38157900E-01, 1.03124990E-01, DVDQ0313
     2  8.14144780E-02, 6.67641930E-02, 5.62453730E-02, 4.83369670E-02, DVDQ0314
     3  4.21809010E-02, 3.72591170E-02, 3.32406510E-02, 2.99038080E-02, DVDQ0315
     4  2.70941180E-02, 2.47004450E-02, 2.26406010E-02, 2.08524830E-02/ DVDQ0316
      DATA (ETA(I, 3), I=1,18)/  1.42857140E-01, 2.85714280E-01,        DVDQ0317
     1  3.42857140E-01, 3.46153840E-01, 2.45614040E-01, 1.87500000E-01, DVDQ0318
     2  1.50303650E-01, 1.24626490E-01, 1.05873640E-01, 9.15858320E-02, DVDQ0319
     3  8.03445710E-02, 7.12783130E-02, 6.38220510E-02, 5.75925170E-02, DVDQ0320
     4  5.23196800E-02, 4.78073100E-02, 4.39090440E-02, 4.05133970E-02/ DVDQ0321
      DATA (ETA(I, 4), I=1,18)/  1.11111110E-01, 2.22222220E-01,        DVDQ0322
     1  2.85714280E-01, 2.53968250E-01, 3.07017540E-01, 2.50000000E-01, DVDQ0323
     2  2.08755060E-01, 1.78037850E-01, 1.54399060E-01, 1.35682710E-01, DVDQ0324
     3  1.20516850E-01, 1.07997450E-01, 9.75059080E-02, 8.86038720E-02, DVDQ0325
     4  8.09709320E-02, 7.43669270E-02, 6.86078780E-02, 6.35504300E-02/ DVDQ0326
      DATA (ETA(I, 5), I=1,18)/  9.09090910E-02, 1.81818180E-01,        DVDQ0327
     1  2.42424240E-01, 2.42424240E-01, 1.73160170E-01, 2.50000000E-01, DVDQ0328
     2  2.27732800E-01, 2.05428290E-01, 1.85278880E-01, 1.67608050E-01, DVDQ0329
     3  1.52231820E-01, 1.38853860E-01, 1.27181620E-01, 1.16957100E-01, DVDQ0330
     4  1.07961240E-01, 1.00010690E-01, 9.29526080E-02, 8.66596710E-02/ DVDQ0331
      DATA (ETA(I, 6), I=1,18)/  7.69230760E-02, 1.53846150E-01,        DVDQ0332
     1  2.09790210E-01, 2.23776220E-01, 1.86480190E-01, 1.11888110E-01, DVDQ0333
     2  1.91295550E-01, 1.91733070E-01, 1.85278880E-01, 1.75988460E-01, DVDQ0334
     3  1.65763530E-01, 1.55516330E-01, 1.45680770E-01, 1.36449950E-01, DVDQ0335
     4  1.27892540E-01, 1.20012830E-01, 1.12782500E-01, 1.06158090E-01/ DVDQ0336
      DATA (ETA(I, 7), I=1,18)/  6.66666660E-02, 1.33333330E-01,        DVDQ0337
     1  1.84615380E-01, 2.05128200E-01, 1.86480190E-01, 1.34265730E-01, DVDQ0338
     2  6.96192690E-02, 1.39442230E-01, 1.52023690E-01, 1.56434190E-01, DVDQ0339
     3  1.56012730E-01, 1.52787970E-01, 1.47993160E-01, 1.42382550E-01, DVDQ0340
     4  1.36418720E-01, 1.30384310E-01, 1.24449660E-01, 1.18714420E-01/ DVDQ0341
      DATA (ETA(I, 8), I=1,18)/  5.88235290E-02, 1.17647060E-01,        DVDQ0342
     1  1.64705880E-01, 1.88235290E-01, 1.80995470E-01, 1.44796380E-01, DVDQ0343
     2  9.21431500E-02, 4.21225830E-02, 9.77295190E-02, 1.14931240E-01, DVDQ0344
     3  1.25367380E-01, 1.30961120E-01, 1.33193850E-01, 1.33136920E-01, DVDQ0345
     4  1.31546610E-01, 1.28951510E-01, 1.25719550E-01, 1.22106250E-01/ DVDQ0346
      DATA (ETA(I, 9), I=1,18)/  5.26315790E-02, 1.05263160E-01,        DVDQ0347
     1  1.48606810E-01, 1.73374610E-01, 1.73374610E-01, 1.48606810E-01, DVDQ0348
     2  1.06692070E-01, 6.09668970E-02, 2.49410030E-02, 6.63064840E-02, DVDQ0349
     3  8.35782530E-02, 9.62949390E-02, 1.05153040E-01, 1.10947430E-01, DVDQ0350
     4  1.14388350E-01, 1.16056360E-01, 1.16406990E-01, 1.15790400E-01/ DVDQ0351
      DATA (ETA(I,10), I=1,18)/  4.76190480E-02, 9.52380950E-02,        DVDQ0352
     1  1.35338350E-01, 1.60401000E-01, 1.65118680E-01, 1.48606810E-01, DVDQ0353
     2  1.15583080E-01, 7.54828240E-02, 3.91930050E-02, 1.45159280E-02, DVDQ0354
     3  4.37790860E-02, 5.88469080E-02, 7.14002090E-02, 8.13614340E-02, DVDQ0355
     4  8.89687070E-02, 9.45644350E-02, 9.84982250E-02, 1.01086850E-01/ DVDQ0356
      DATA (ETA(I,11), I=1,18)/  4.34782610E-02, 8.69565210E-02,        DVDQ0357
     1  1.24223600E-01, 1.49068320E-01, 1.56914020E-01, 1.46453090E-01, DVDQ0358
     2  1.20608430E-01, 8.61488760E-02, 5.16893250E-02, 2.46139640E-02, DVDQ0359
     3  8.33088030E-03, 2.82465160E-02, 4.03201180E-02, 5.13861540E-02, DVDQ0360
     4  6.10071000E-02, 6.90731460E-02, 7.56466410E-02, 8.08694600E-02/ DVDQ0361
      DATA (ETA(I,12), I=1,18)/  4.00000000E-02, 7.99999990E-02,        DVDQ0362
     1  1.14782610E-01, 1.39130430E-01, 1.49068320E-01, 1.43105590E-01, DVDQ0363
     2  1.23020590E-01, 9.37299770E-02, 6.20271900E-02, 3.44595500E-02, DVDQ0364
     3  1.51622020E-02, 4.72588120E-03, 1.78691420E-02, 2.69906870E-02, DVDQ0365
     4  3.60496320E-02, 4.45264850E-02, 5.21503400E-02, 5.88141320E-02/ DVDQ0366
      DATA (ETA(I,13), I=1,18)/  3.70370370E-02, 7.40740740E-02,        DVDQ0367
     1  1.06666670E-01, 1.30370370E-01, 1.41706920E-01, 1.39130430E-01, DVDQ0368
     2  1.23671500E-01, 9.89371980E-02, 7.02974820E-02, 4.33935080E-02, DVDQ0369
     3  2.24625220E-02, 9.18921340E-03, 2.65466170E-03, 1.11137880E-02, DVDQ0370
     4  1.77085690E-02, 2.47369290E-02, 3.17436890E-02, 3.84252050E-02/ DVDQ0371
      DATA (ETA(I,14), I=1,18)/  3.44827580E-02, 6.89655170E-02,        DVDQ0372
     1  9.96168580E-02, 1.22605360E-01, 1.34865900E-01, 1.34865900E-01, DVDQ0373
     2  1.23138430E-01, 1.02348820E-01, 7.67616190E-02, 5.11744120E-02, DVDQ0374
     3  2.96272910E-02, 1.43647470E-02, 5.49240340E-03, 1.47872400E-03, DVDQ0375
     4  6.81096220E-03, 1.14170350E-02, 1.66487940E-02, 2.21683550E-02/ DVDQ0376
      DATA (ETA(I,15), I=1,18)/  3.22580640E-02, 6.45161290E-02,        DVDQ0377
     1  9.34371530E-02, 1.15684090E-01, 1.28537880E-01, 1.30515390E-01, DVDQ0378
     2  1.21814360E-01, 1.04412310E-01, 8.17139820E-02, 5.77775630E-02, DVDQ0379
     3  3.63173250E-02, 1.98094500E-02, 9.03588950E-03, 3.24365270E-03, DVDQ0380
     4  8.17727570E-04, 4.12042260E-03, 7.24846120E-03, 1.10152930E-02/ DVDQ0381
      DATA (ETA(I,16), I=1,18)/  3.03030310E-02, 6.06060610E-02,        DVDQ0382
     1  8.79765400E-02, 1.09481920E-01, 1.22695250E-01, 1.26200830E-01, DVDQ0383
     2  1.19968690E-01, 1.05466980E-01, 8.54282550E-02, 6.32801880E-02, DVDQ0384
     3  4.23702130E-02, 2.52120280E-02, 1.30062050E-02, 5.60267280E-03, DVDQ0385
     4  1.89564120E-03, 4.49337170E-04, 2.46448470E-03, 4.53959090E-03/ DVDQ0386
      DATA (ETA(I,17), I=1,18)/  2.85714290E-02, 5.71428580E-02,        DVDQ0387
     1  8.31168840E-02, 1.03896110E-01, 1.17302050E-01, 1.21994140E-01, DVDQ0388
     2  1.17787440E-01, 1.05768320E-01, 8.81402630E-02, 6.78002020E-02, DVDQ0389
     3  4.77313430E-02, 3.03744910E-02, 1.71681900E-02, 8.40400930E-03, DVDQ0390
     4  3.43020790E-03, 1.09766650E-03, 2.45530670E-04, 1.45910240E-03/ DVDQ0391
      DATA (ETA(I,18), I=1,18)/  2.70270270E-02, 5.40540540E-02,        DVDQ0392
     1  7.87644800E-02, 9.88417010E-02, 1.12320110E-01, 1.17936120E-01, DVDQ0393
     2  1.15399860E-01, 1.05508440E-01, 9.00459990E-02, 7.14650780E-02, DVDQ0394
     3  5.24077240E-02, 3.51828080E-02, 2.13442370E-02, 1.14930500E-02, DVDQ0395
     4  5.35390550E-03, 2.07666640E-03, 6.30416580E-04, 1.33499980E-04/ DVDQ0396
C     ETA(I,J) I=1,2,...,J IS USED IN THE FIRST MODIFICATION OF THE     DVDQ0397
C     I-TH DIFFERENCE OF A J-TH ORDER METHOD AFTER THE STEPSIZE IS      DVDQ0398
C     HALVED.                                                           DVDQ0399
C     ETA(I,J) J=1,2,...,I-1 IS USED IN THE SECOND MODIFICATION OF      DVDQ0400
C     THE (J+1)-ST DIFFERENCE OF AN I-TH ORDER METHOD.                  DVDQ0401
C     THE TWO MODIFICATIONS OF THE DIFFERENCE TABLE AFTER HALVING THE   DVDQ0402
C     STEPSIZE REMOVES MOST OF THE INSTABILITY INHERENT IN THE METHOD   DVDQ0403
C     USED HERE FOR HALVING THE STEPSIZE.                               DVDQ0404
C                                                                       DVDQ0405
C                                                                       DVDQ0406
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE FOLLOWING CARD.    DVDQ0407
      DATA LGSS,LGSD,LGSE/0,0,0/                                        DVDQ0408
C                                                                       DVDQ0409
C****                                                                           
C     THE STATEMENTS BETWEEN THE LINE ABOVE AND A SIMILAR LINE BELOW            
C     ARE INSERTED IN THE JPL VERSION OF THIS SUBROUTINE TO INSURE              
C     THAT THE CODE IS COMPILED CORRECTLY. THESE STATEMENTS THEMSELVES          
C     ARE NOT ACTUALLY COMPILED.                                                
C                                                                               
C      NEQ=NEQ                                                                   
C      HMINA=HMINA                                                               
C      HMAXA=HMAXA                                                               
C      DELT=DELT                                                                 
C      TFINAL=TFINAL                                                             
C      MXSTEP=MXSTEP                                                             
C****                                                                           
C                                                                       DVDQ0410
      SAVE
C     INITIALIZE                                                        DVDQ0411
      KSTEP=-1                                                          DVDQ0412
      NE=NEQ                                                            DVDQ0413
      IF (NE.LE.0) GO TO 1190                                           DVDQ0414
      HH=H                                                              DVDQ0415
      NV=0                                                              DVDQ0416
      KDMAX=0                                                           DVDQ0417
      KDD=KD(1)                                                         DVDQ0418
      KDS=KDD                                                           DVDQ0419
      DO 10 I=1,NE                                                      DVDQ0420
      KQ(I)=1                                                           DVDQ0421
      DT(1,I)=0.D0                                                      DVDQ0422
      IF (KDS.LE.0) KDD=IABS(KD(I))                                     DVDQ0423
      IF ((KDD.EQ.0).OR.(KDD.GT.KMAXO)) HH=0.D0                         DVDQ0424
      IF (KDD.GT.KDMAX) KDMAX=KDD                                       DVDQ0425
   10 NV=NV+KDD                                                         DVDQ0426
C                                                                       DVDQ0427
      IF ((DELT*HH).LE.0.D0) GO TO 1190                                 DVDQ0428
      ERRMX=P1                                                          DVDQ0429
      EMAX=ERND                                                         DVDQ0430
      RNDC=RND*P25                                                      DVDQ0431
      LDOUB=0                                                           DVDQ0432
      E2HFAC=P25                                                        DVDQ0433
      LSC=8                                                             DVDQ0434
      LSTC=4                                                            DVDQ0435
C     LSC AND LSTC ARE USED IN COMBINATION AS FOLLOWS                   DVDQ0436
C       LSTC=4,  LSC=4     FIRST TIME THROUGH THE FIRST STEP            DVDQ0437
C       LSTC=3,  LSC=4     SECOND TIME THROUGH THE FIRST STEP           DVDQ0438
C                          (NECESSARY TO CHECK STABILITY)               DVDQ0439
C       LSTC=2,  LSC=4     THIRD TIME THROUGH THE FIRST STEP            DVDQ0440
C                          (ONLY OCCURS IF INSTABILITY POSSIBLE)        DVDQ0441
C       LSTC=2,  LSC=2     SECOND STEP  (IF KQ(I)=2 , I=1,...,NEQ)      DVDQ0442
C       LSTC=1,  LSC=0     STARTING, ONE DERIVATIVE EVAL. PER STEP.     DVDQ0443
C       LSTC=1,  LSC.GT.0  SET WHEN STARTING TWO DERIV. EVAL. PER STEP  DVDQ0444
C       LSTC=-1  LSC.LT.0  SET WHEN HALVING THE STEPSIZE                DVDQ0445
C      IN THE LAST TWO CASES LSC IS SET EQUAL TO LSTC*(MAXIMUM KQ(I)    DVDQ0446
C      +1). AT THE END OF EACH STEP IF LSC.NE.0 IT IS REPLACED BY       DVDQ0447
C      LSC-LSTC UNTIL LSC=0, AT WHICH TIME LSTC IS SET EQUAL TO 0.      DVDQ0448
C       WHEN DOUBLING H, LSTC IS SET EQUAL TO -1 AND LSC TO -3.         DVDQ0449
C      UNDER CERTAIN CONDITIONS WHEN KQ(I)=1, LSTC IS SET =-1 AND LSC=-5DVDQ0450
C                                                                       DVDQ0451
      KSOUT=MXSTEP                                                      DVDQ0452
      TOUT=T                                                            DVDQ0453
      IFL=13                                                            DVDQ0454
   20 IFLAG=1                                                           DVDQ0455
      GO TO 315                                                         DVDQ0456
C     END OF INITIALIZATION                                             DVDQ0457
C                                                                       DVDQ0458
C                                                                       DVDQ0459
      ENTRY DVDQ1 (NEQ,KD,IFLAG,MXSTEP,KSTEP,KEMAX,KQ,HMINA,HMAXA,      DVDQ0460
     1EMAX,T,Y,F,H,DELT,TFINAL,YN,DT,EP)                                DVDQ0461
C                                                                       DVDQ0462
C     TO OUTPUT VARIABLES IN THE CALLING SEQUENCE REMOVE THE C-S        DVDQ0463
C     IN COLUMN ONE OF THE FOLLOWING CARDS UNTIL REACHING THE COMMENT   DVDQ0464
C     END OF CODE FOR PRINTING VARIABLES IN CALLING SEQUENCE.           DVDQ0465
C     IF (NEQ.NE.0) GO TO 28                                            DVDQ0466
C     NEQ=1                                                             DVDQ0467
C  22 WRITE(6,5000) T,DELT,HMINA,HMAXA,KEMAX,EMAX,IFLAG,TFINAL,MXSTEP   DVDQ0468
C5000 FORMAT(3H0T=1PD24.17,7H  DELT=D12.5,8H  HMINA=,E10.3,8H  HMAXA=,  DVDQ0469
C    1  E10.3,8H  KEMAX=,I2,7H  EMAX=E10.3,8H  IFLAG=,I2/               DVDQ0470
C    2  9H  I KQ KD,7X,4HF(I),9X,1HJ,12X,4HY(J),22X,5HYN(J),            DVDQ0471
C    3  10X,7HTFINAL=1PD15.8,9H  MXSTEP=I4)                             DVDQ0472
C     J=1                                                               DVDQ0473
C     DO 24 I=1,NE                                                      DVDQ0474
C     IF (KDS.LT.0) KDD=IABS(KD(I))                                     DVDQ0475
C     K=KDD                                                             DVDQ0476
C     WRITE(6,5001) I,KQ(I),KDD,F(I),J,Y(J),YN(J)                       DVDQ0477
C5001 FORMAT(1H ,I2,2I3,1PD17.8,I4,2D26.17)                             DVDQ0478
C  23 J=J+1                                                             DVDQ0479
C     K=K-1                                                             DVDQ0480
C     IF (K.EQ.0) GO TO 24                                              DVDQ0481
C     WRITE(6,5002) J,Y(J),YN(J)                                        DVDQ0482
C5002 FORMAT(26X,I4,1P2D26.17)                                          DVDQ0483
C     GO TO 23                                                          DVDQ0484
C  24 CONTINUE                                                          DVDQ0485
C     WRITE(6,5003)                                                     DVDQ0486
C5003 FORMAT(3H0 I,15X,16HDIFFERENCE TABLE)                             DVDQ0487
C     DO 27 I=1,NE                                                      DVDQ0488
C     KQQ=KQ(I)+1                                                       DVDQ0489
C     K=MIN0(KQQ,7)                                                     DVDQ0490
C     WRITE(6,5004) I,(DT(IO,I),IO=1,K)                                 DVDQ0491
C5004 FORMAT(1H ,I2,1PD19.8,6D16.7)                                     DVDQ0492
C                                                                       DVDQ0493
C     IF (K.EQ.KQQ) GO TO 27                                            DVDQ0494
C     K=K+1                                                             DVDQ0495
C     WRITE(6,5005) (DT(IO,I),IO=K,KQQ)                                 DVDQ0496
C5005 FORMAT(1H ,1PD21.5,7D14.5)                                        DVDQ0497
C  27 CONTINUE                                                          DVDQ0498
C     IF (NEQ.EQ.0) RETURN                                              DVDQ0499
C     NEQ=0                                                             DVDQ0500
C 28  CONTINUE                                                          DVDQ0501
C     END OF CODE FOR PRINTING VARIABLES IN CALLING SEQUENCE.           DVDQ0502
C                                                                       DVDQ0503
      IF (2-IFL) 30,60,320                                              DVDQ0504
   30 IF (IFL.GT.5) GO TO 1180                                          DVDQ0505
C                                                                       DVDQ0506
C     SET STEP STOP                                                     DVDQ0507
      KSOUT=KSTEP+MXSTEP                                                DVDQ0508
      IF (IFL-4) 40,1210,210                                            DVDQ0509
C                                                                       DVDQ0510
C                                                                       DVDQ0511
C     SET PRINT STOP                                                    DVDQ0512
   40 TOUT=T+DELT                                                       DVDQ0513
C                                                                       DVDQ0514
   50 TPS1=ABS(SNGL(DMOD((TOUT-TL)/HH,2.D0))-PTS1)                      DVDQ0515
      LFD=-1                                                            DVDQ0516
      IF (TPS1.GE.P5) LFD=1                                             DVDQ0517
C                                                                       DVDQ0518
C     LFD IS USED TO INDICATE WHETHER DOUBLING H IS PERMITTED.          DVDQ0519
C     IF LFD.LT.0 AT THE END OF A STEP THEN DOUBLING H IS               DVDQ0520
C     NOT PERMITTED. THE SIGN OF LFD IS CHANGED JUST BEFORE THE         DVDQ0521
C     END OF EACH STEP. IF DELT=H*(POWER OF 2) THEN                     DVDQ0522
C     OUTPUT VALUES WILL BE OBTAINED WITHOUT INTERPOLATION.             DVDQ0523
C                                                                       DVDQ0524
      GO TO 200                                                         DVDQ0525
C                                                                       DVDQ0526
C                                                                       DVDQ0527
C     ENTRY WITH IFLAG=2                                                DVDQ0528
C                                                                       DVDQ0529
C     UPDATE DIFFERENCE TABLE                                           DVDQ0530
C     AND COMPUTE KQM=MAXIMUM VALUE OF KQ(I), I=1,2,...,NEQ.            DVDQ0531
C                                                                       DVDQ0532
   60 KQM=0                                                             DVDQ0533
      DO 80 I=1,NE                                                      DVDQ0534
      KQQ=KQ(I)                                                         DVDQ0535
      IF (KQQ.GT.KQM) KQM=KQQ                                           DVDQ0536
      D(1)=F(I)                                                         DVDQ0537
      DO 70 K=1,KQQ                                                     DVDQ0538
      D(K+1)=D(K)-DT(K,I)                                               DVDQ0539
   70 DT(K,I)=D(K)                                                      DVDQ0540
      DT(KQQ+1,I)=D(KQQ+1)                                              DVDQ0541
   80 CONTINUE                                                          DVDQ0542
C     END OF UPDATING DIFFERENCE TABLE                                  DVDQ0543
C                                                                       DVDQ0544
C     STORE Y(J) IN YN(J)                                               DVDQ0545
      DO 90 J=1,NV                                                      DVDQ0546
   90 YN(J)=Y(J)                                                        DVDQ0547
C                                                                       DVDQ0548
      LFD=-LFD                                                          DVDQ0549
      TL=T                                                              DVDQ0550
      KSTEP=KSTEP+1                                                     DVDQ0551
C                                                                       DVDQ0552
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE 2 FOLLOWING CARDS. DVDQ0553
      IF (LGSS) 1430,110,1510                                           DVDQ0554
  100 IFLAG=2                                                           DVDQ0555
  110 IF (LSC.EQ.0) GO TO 140                                           DVDQ0556
      LSC=LSC-LSTC                                                      DVDQ0557
      IF (LSC.EQ.0) GO TO 130                                           DVDQ0558
      IF (LSTC.NE.(-1)) GO TO 140                                       DVDQ0559
      IF (LDOUB.LT.0) RNDC=RND*P1                                       DVDQ0560
  120 E2HAVE=E2HMAX                                                     DVDQ0561
      TPS1=PTS1                                                         DVDQ0562
      GO TO 190                                                         DVDQ0563
  130 IF (ABS(SNGL(HH)).LT.HMINA) GO TO 1000                            DVDQ0564
      LSTC=0                                                            DVDQ0565
  140 IF (LDOUB.NE.1) GO TO 150                                         DVDQ0566
      IF ((LFD.GT.0).AND.(ABS(SNGL(HH+HH)).LE.HMAXA)) GO TO 1030        DVDQ0567
      GO TO 200                                                         DVDQ0568
  150 RQMAX=PTS1/FLOAT(KQM+3)                                           DVDQ0569
      IF (LSTC.NE.0.OR.E2HAVE.EQ.0.E0)       GO TO 120                        FK
      TPS1=E2HMAX/E2HAVE                                                DVDQ0571
      IF (TPS1-PTS1) 160,190,170                                        DVDQ0572
  160 E2HFAC=AMAX1(.075E0,E2HFAC-RQMAX,E2HFAC*TPS1)                     DVDQ0573
      GO TO 180                                                         DVDQ0574
  170 TPS1=TPS1*TPS1                                                    DVDQ0575
      E2HFAC=AMIN1(PTS1,E2HFAC*TPS1)                                    DVDQ0576
  180 RNDC=(1.1-E2HFAC)*RND                                             DVDQ0577
      E2HAVE=P5*(E2HMAX+E2HAVE)                                         DVDQ0578
  190 ERRMX=AMAX1(P1,ERRMX-RQMAX*TPS1)                                  DVDQ0579
C     E2HFAC IS A FACTOR WHICH IS TAKEN TIMES AN INITIAL ESTIMATE OF    DVDQ0580
C            E2H TO GET A FINAL VALUE OF E2H. (E2H=ESTIMATE OF WHAT     DVDQ0581
C            (ESTIMATED ERROR)/(REQUESTED ERROR) WOULD BE IF H WERE     DVDQ0582
C            DOUBLED.)                                                  DVDQ0583
C     E2HMAX IS THE MAXIMUM VALUE OF THE INITIAL ESTIMATE OF E2H OVER   DVDQ0584
C            ALL COMPONENTS WITH KQ(I).GT.1.                            DVDQ0585
C     E2HAVE IS A WEIGHTED AVERAGE OF PAST VALUES OF E2HMAX.            DVDQ0586
C     THE VALUE OF E2HFAC TENDS TO BE SMALLER WHEN E2HMAX IS            DVDQ0587
C     CONSISTANTLY SMALLER THAN E2HAVE.                                 DVDQ0588
C                                                                       DVDQ0589
C                                                                       DVDQ0590
C     CHECK FOR PRINT STOP AND FOR T REACHING TFINAL                    DVDQ0591
  200 TPD=(TOUT-TL)/HH                                                  DVDQ0592
      TPD1=(TFINAL-TL)/HH                                               DVDQ0593
C                                                                       DVDQ0594
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE FOLLOWING CARD.    DVDQ0595
      IF (LGSE.LT.0) GO TO 1780                                         DVDQ0596
      IF (TPD1.LT.FAC(1)) GO TO 1220                                    DVDQ0597
      IF (TPD.LE.0.D0) GO TO 1280                                       DVDQ0598
C                                                                       DVDQ0599
C     CHECK FOR STEP STOP                                               DVDQ0600
      IF (KSOUT.GT.KSTEP) GO TO 210                                     DVDQ0601
C                                                                       DVDQ0602
      IFL=5                                                             DVDQ0603
      GO TO 310                                                         DVDQ0604
C                                                                       DVDQ0605
C     CHECK TO SEE IF ROUND-OFF ERROR IS PROMINENT                      DVDQ0606
  210 IF (EMAX.EQ.ERND) GO TO 220                                       DVDQ0607
C     IT IS                                                             DVDQ0608
      IFL=6                                                             DVDQ0609
      IF (EMAX.GE.P1) GO TO 310                                         DVDQ0610
      IF ((LSTC.GE.0).OR.(LDOUB.EQ.1)) ERRMX=PTS1                       DVDQ0611
C                                                                       DVDQ0612
  220 IFL=1                                                             DVDQ0613
  230 T=TL+HH                                                           DVDQ0614
C                                                                       DVDQ0615
C     START A NEW STEP                                                  DVDQ0616
C                                                                       DVDQ0617
C     PREDICT                                                           DVDQ0618
  240 J=0                                                               DVDQ0619
      DO 290 I=1,NE                                                     DVDQ0620
      IF (KDS.LE.0) KDD=IABS(KD(I))                                     DVDQ0621
      KDC=KDD                                                           DVDQ0622
  250 KQQ=KQ(I)                                                         DVDQ0623
      TPD=0.D0                                                          DVDQ0624
      K=KDC                                                             DVDQ0625
  260 TPD=TPD+DT(KQQ,I)*GAM(KQQ,KDC)                                    DVDQ0626
      KQQ=KQQ-1                                                         DVDQ0627
      IF (KQQ.GT.0) GO TO 260                                           DVDQ0628
  270 K=K-1                                                             DVDQ0629
      IF (K.LE.0) GO TO 280                                             DVDQ0630
      L=J+K                                                             DVDQ0631
      TPD=YN(L+1)*FAC(K)+HH*TPD                                         DVDQ0632
      GO TO 270                                                         DVDQ0633
  280 J=J+1                                                             DVDQ0634
      Y(J)=YN(J)+HH*TPD                                                 DVDQ0635
      KDC=KDC-1                                                         DVDQ0636
      IF (KDC.GT.0) GO TO 250                                           DVDQ0637
  290 CONTINUE                                                          DVDQ0638
C     END OF PREDICT                                                    DVDQ0639
C                                                                       DVDQ0640
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C  IN COLUMN ONE   DVDQ0641
C     OF THE 2 FOLLOWING CARDS                                          DVDQ0642
C     IF (IFL) 20,320,300                                               DVDQ0643
C 300 CONTINUE                                                          DVDQ0644
C     AND THEN REMOVE THE 2 FOLLOWING CARDS.                            DVDQ0645
      IF (IFL) 1240,320,300                                             DVDQ0646
  300 IF (LGSD.NE.0) GO TO 1520                                         DVDQ0647
C                                                                       DVDQ0648
  310 IFLAG=IFL                                                         DVDQ0649
  315 CONTINUE                                                          DVDQ0650
C                                                                       DVDQ0651
C     TO OUTPUT VARIABLES IN THE CALLING SEQUENCE REMOVE THE C IN       DVDQ0652
C     COLUMN ONE OF THE FOLLOWING CARD.                                 DVDQ0653
C     IF (NEQ.EQ.0) GO TO 22                                            DVDQ0654
C                                                                       DVDQ0655
      RETURN                                                            DVDQ0656
C                                                                       DVDQ0657
C                                                                       DVDQ0658
C     ENTRY WITH IFLAG=1                                                DVDQ0659
  320 EPS=EP(1)                                                         DVDQ0660
      ERND=0.                                                           DVDQ0661
      EMAX=0.                                                           DVDQ0662
      E2HMAX=0.                                                         DVDQ0663
      J=0                                                               DVDQ0664
      IF (LDOUB.GE.0) LDOUB=1                                           DVDQ0665
C                                                                       DVDQ0666
C     LDOUB IS SET IN THE LOOP BELOW AS FOLLOWS                         DVDQ0667
C     LDOUB=0   HALVE                                                   DVDQ0668
C     LDOUB=1   DOUBLE                                                  DVDQ0669
C     LDOUB=2   DO NOT DOUBLE                                           DVDQ0670
C                                                                       DVDQ0671
C     LDOUB,LT.0 AT THE BEGINNING OF THE LOOP INDICATES THE FOLLOWING   DVDQ0672
C       =-3  STEPSIZE HAS JUST BEEN HALVED. IF A DISCONTINUITY IS       DVDQ0673
C            NOT INDICATED MODIFY THE DIFFERENCE TABLE AND REPEAT       DVDQ0674
C            THE STEP.                                                  DVDQ0675
C       =-2  STEP AFTER LDOUB=-3. PROCEED AS USUAL (ORDER IS NOT        DVDQ0676
C            CHANGED)                                                   DVDQ0677
C       =-1  STEP AFTER LDOUB=-2. MODIFY THE DIFFERENCE TABLE ONCE      DVDQ0678
C            AGAIN AND REPEAT THE STEP.                                 DVDQ0679
C     IF LDOUB IS SET EQUAL TO -4 THE ORDER IN AT LEAST ONE COMPONENT   DVDQ0680
C     HAS BEEN GREATLY REDUCED AND THE STEP IS REPEATED.                DVDQ0681
C                                                                       DVDQ0682
C                                                                       DVDQ0683
C     IF THE OUTPUT OPTION IS ELIMINATED, REMOVE THE 4 FOLLOWING CARDS. DVDQ0684
      IF (NEQ.LE.0) WRITE (6,5020) LSC,LFD,LSTC,KSTEP,E2HFAC,ERRMX,HH   DVDQ0685
 5020 FORMAT (19H0 I  KQQ LRND LDOUB,5X,1HE,9X,3HE2H,                   DVDQ0686
     1  8X,3HEPS,3X,4HLSC=,I3,6H  LFD=,I2,7H  LSTC=,I2,8H  KSTEP=,I4,   DVDQ0687
     2  9H  E2HFAC=,F4.2,8H  ERRMX=,F4.2,4H  H=,1PD9.2)                 DVDQ0688
C                                                                       DVDQ0689
C                                                                       DVDQ0690
C     BEGINNING OF LOOP FOR CORRECTING, ESTIMATING THE ERROR,           DVDQ0691
C     AND ADJUSTING THE NUMBER OF DIFFERENCES USED                      DVDQ0692
C                                                                       DVDQ0693
      DO 790 I=1,NE                                                     DVDQ0694
      IF (KDS.LE.0) KDD=IABS(KD(I))                                     DVDQ0695
      KQQ=KQ(I)                                                         DVDQ0696
C     KQQ GIVES THE ORDER OF THE PREDICTOR FORMULA AND KQQ+1 THE        DVDQ0697
C     ORDER OF THE CORRECTOR FORMULA.                                   DVDQ0698
C                                                                       DVDQ0699
      KQ1=KQQ+1                                                         DVDQ0700
      D(1)=F(I)                                                         DVDQ0701
C     FORM THE DIFFERENCE TABLE FROM PREDICTED DERIVATIVE VALUES.       DVDQ0702
      DO 330 K=1,KQ1                                                    DVDQ0703
      D(K+1)=D(K)-DT(K,I)                                               DVDQ0704
  330 CONTINUE                                                          DVDQ0705
C     D(K) GIVES THE (K-1)-ST DIFFERENCE FORMED FROM PREDICTED          DVDQ0706
C     DERIVATIVE VALUES                                                 DVDQ0707
      TPS3=ABS(SNGL(D(KQQ+1)))                                          DVDQ0708
      IF (LDOUB.LT.0) GO TO 720                                         DVDQ0709
C                                                                       DVDQ0710
  340 IF (KQQ.NE.1) GO TO 520                                           DVDQ0711
C                                                                       DVDQ0712
C     KQ(I)=1  IS TREATED AS A SPECIAL CASE                             DVDQ0713
      E2H=PTS2                                                          DVDQ0714
      TPS5=DT(3,I)                                                      DVDQ0715
      IF (LSTC.LT.2) GO TO 370                                          DVDQ0716
C     FIRST STEP OF INTEGRATION                                         DVDQ0717
      IF (LSTC.NE.4) GO TO 350                                          DVDQ0718
      TPS4=0.                                                           DVDQ0719
      IF (KDD.GT.1) TPS3=AMAX1(TPS3,ABS(SNGL(HH*D(1))))                 DVDQ0720
      TPS3=TPS3*P1                                                      DVDQ0721
      GO TO 510                                                         DVDQ0722
  350 DT(2,I)=D(2)                                                      DVDQ0723
      D(2)=D(1)-DT(5,I)                                                 DVDQ0724
      TPS2=-D(2)                                                        DVDQ0725
      TPS3=PTS5*ABS(TPS2)                                               DVDQ0726
C     FIRST STEP THAT KQ(I)=1                                           DVDQ0727
  360 DT(7,I)=PT(4)                                                     DVDQ0728
      IF (LSTC-2) 420,380,380                                           DVDQ0729
  370 IF (TPS5.EQ.0.) GO TO 360                                         DVDQ0730
      IF (DT(6,I).EQ.0.D0) GO TO 400                                    DVDQ0731
      TPS2=DT(5,I)-DT(1,I)                                              DVDQ0732
  380 TPS4=DT(4,I)                                                      DVDQ0733
      TPS1=ABS(TPS4)                                                    DVDQ0734
      TPS4=TPS2*SIGN(PTS2,TPS4)-TPS5*TPS1                               DVDQ0735
      IF (TPS4.GT.(-TPS1)) GO TO 410                                    DVDQ0736
  390 TPS6=-PTS1                                                        DVDQ0737
      GO TO 450                                                         DVDQ0738
C     FIRST STEP AFTER THE STEPSIZE HAS BEEN CHANGED                    DVDQ0739
  400 DT(6,I)=PT(1)                                                     DVDQ0740
      TPS6=0.                                                           DVDQ0741
      GO TO 450                                                         DVDQ0742
  410 IF (TPS4.LT.TPS1) GO TO 440                                       DVDQ0743
      IF (TPS1.EQ.0.) GO TO 390                                         DVDQ0744
  420 TPS6=PTS1                                                         DVDQ0745
      GO TO 450                                                         DVDQ0746
  430 KQ(I)=2                                                           DVDQ0747
      IF (2-LSTC) 510,510,520                                           DVDQ0748
  440 TPS6=TPS4/TPS1                                                    DVDQ0749
  450 TPS4=TPS5+TPS6                                                    DVDQ0750
      IF (TPS4.LT.P25) GO TO 430                                        DVDQ0751
C     INCREASE E2H IF (-S).GT..25                                       DVDQ0752
      E2H=PTS4*TPS4                                                     DVDQ0753
      IF (2-LSTC) 460,470,480                                           DVDQ0754
  460 LSC=0                                                             DVDQ0755
      GO TO 510                                                         DVDQ0756
  470 IF (TPS5-P25) 430,460,460                                         DVDQ0757
  480 IF (TPS4.GT.PTS2) GO TO 490                                       DVDQ0758
      IF (TPS4.GT.P5) D(2)=D(2)*GAM(2,1)                                DVDQ0759
      GO TO 510                                                         DVDQ0760
  490 IF (TPS4.LT.PTS4) GO TO 500                                       DVDQ0761
      TPS4=PTS4                                                         DVDQ0762
      D(2)=D(2)/PT(3)                                                   DVDQ0763
C     THE ESTIMATE OF E (AND HENCE OF E2H) IS INCREASED IF (-S).GE.8.   DVDQ0764
      TPS3=TPS3*SNGL(DT(7,I))                                           DVDQ0765
      GO TO 510                                                         DVDQ0766
  500 D(2)=D(2)*DBLE(PTS2*(TPS4-PTS1)/(TPS4*TPS4))                      DVDQ0767
      IF (TPS4.GE.P3E1) E2H=E2H*SNGL(DT(7,I))                           DVDQ0768
C     STORE D(1)=PREDICTED DERIVATIVE  AND  D(2)=2*(CORRECTED Y -       DVDQ0769
C     PREDICTED Y)/H    D(1) AND D(2) ARE USED TO COMPUTE (-S) ON       DVDQ0770
C     THE NEXT STEP.                                                    DVDQ0771
  510 DT(5,I)=D(1)                                                      DVDQ0772
      DT(4,I)=D(2)                                                      DVDQ0773
      D(4)=TPS4                                                         DVDQ0774
C     STORE D(4)= CURRENT ESTIMATE OF (-S).  (-S).GT.3 IS AN INDICATION DVDQ0775
C     THAT THE STEPSIZE SHOULD BE LIMITED BECAUSE OF STABILITY PROBLEMS.DVDQ0776
C     S=H*(ESTIMATE OF EIGENVALUE OF F)=H*(DIFFERENCE BETWEEN PREDICTED DVDQ0777
C     AND CORRECTED DERIVATIVE VALUES)/(DIFFERENCE BETWEEN PREDICTED    DVDQ0778
C     AND CORRECTED INTEGRALS OF THE DERIVATIVE VALUES)                 DVDQ0779
C     THE TREATMENT OF THE CASE KQ(I)=1 COULD BE IMPROVED BY USING A    DVDQ0780
C     SPECIAL METHOD FOR STIFF EQUATIONS WHEN (-S).GT.3  (MAYBE).       DVDQ0781
C     (THE ENTIRE TREATMENT OF THE CASE KQ(I)=1 IS FAR FROM IDEAL.)     DVDQ0782
      DT(3,I)=D(4)                                                      DVDQ0783
C                                                                       DVDQ0784
C     CORRECT                                                           DVDQ0785
  520 KDC=0                                                             DVDQ0786
      TPD=D(KQ1)                                                        DVDQ0787
      J=J+KDD                                                           DVDQ0788
      K=J                                                               DVDQ0789
  530 TPD=HH*TPD                                                        DVDQ0790
      KDC=KDC+1                                                         DVDQ0791
      Y(K)=Y(K)+GAM(KQQ+1,KDC)*TPD                                      DVDQ0792
      K=K-1                                                             DVDQ0793
      IF (KDC.LT.KDD) GO TO 530                                         DVDQ0794
C     END OF CORRECT                                                    DVDQ0795
C                                                                       DVDQ0796
      IF (EPS) 540,550,560                                              DVDQ0797
  540 EPS=EP(I)                                                         DVDQ0798
      IF (EPS.NE.0.) GO TO 560                                          DVDQ0799
  550 IF (HMAXA) 1190,780,1190                                          DVDQ0800
  560 TPS4=ABS(SNGL(D(KQQ+2)))                                          DVDQ0801
      TPS2=ABS(SNGL(D(KQQ)))                                            DVDQ0802
      TPS6=SNGL(HH)/EPS                                                 DVDQ0803
C                                                                       DVDQ0804
      E=ABS(SNGL(GAS(KQQ+1))*TPS3*TPS6)                                 DVDQ0805
C     E GIVES  ABS((ESTIMATED ERROR)/EPS)                               DVDQ0806
C                                                                       DVDQ0807
      LRND=1                                                            DVDQ0808
C                                                                       DVDQ0809
C     LRND= 1  MEANS NO ROUND-OFF ERROR                                 DVDQ0810
C         = 0  MEANS SOME ROUND-OFF ERROR                               DVDQ0811
C         =-1  MEANS EXTREME ROUND-OFF ERROR                            DVDQ0812
C                                                                       DVDQ0813
      FRND=RNDC*ABS(SNGL(PT(KQQ+2)*D(1)))                               DVDQ0814
C     CHECK TO SEE IF ROUND OFF ERROR IS DOMINANT                       DVDQ0815
      IF ((TPS3+TPS4).GT.FRND) GO TO 570                                DVDQ0816
      LRND=0                                                            DVDQ0817
      IF ((PTS4*TPS2).LT.FRND) LRND=-1                                  DVDQ0818
C                                                                       DVDQ0819
  570 IF (E.LE.ERND) GO TO 590                                          DVDQ0820
      IF (E.LE.EMAX) GO TO 580                                          DVDQ0821
      EMAX=E                                                            DVDQ0822
      KEMAX=I                                                           DVDQ0823
  580 IF (LRND.LE.0) GO TO 590                                          DVDQ0824
      ERND=E                                                            DVDQ0825
      IF (ERND.GT.ERRMX) LDOUB=0                                        DVDQ0826
  590 IF (LDOUB.LE.0) GO TO 780                                         DVDQ0827
      TPS1=ABS(SNGL(DD(KQQ)))                                           DVDQ0828
      TPS5=TPS1                                                         DVDQ0829
      IF (KQQ-2) 600,610,620                                            DVDQ0830
  600 E2H=E*E2H                                                         DVDQ0831
      IF (E2H.LT.P01) GO TO 780                                         DVDQ0832
      IF (SNGL(D(4)).LT.P3E1) GO TO 770                                 DVDQ0833
      LSTC=-1                                                           DVDQ0834
      LSC=-5                                                            DVDQ0835
      GO TO 770                                                         DVDQ0836
  610 TPS1=TPS2                                                         DVDQ0837
      IF (LSTC.NE.2) GO TO 620                                          DVDQ0838
      KQ(I)=3                                                           DVDQ0839
      TPS2=0.                                                           DVDQ0840
      TPS4=0.                                                           DVDQ0841
      LRND=0                                                            DVDQ0842
  620 E2H=TPS2+TPS3+TPS4                                                DVDQ0843
      E2H=ABS(SNGL(GAS(KQQ-1)*PT(KQQ+1))*E2H*TPS6)                      DVDQ0844
C     E2H IS USED AS AN ESTIMATE OF WHAT THE VALUE OF E WOULD BE        DVDQ0845
C     IF H WERE DOUBLED. THE ESTIMATE IS CONSERVATIVELY LARGE.          DVDQ0846
      IF (E2H.GT.E2HMAX) E2HMAX=E2H                                     DVDQ0847
C                                                                       DVDQ0848
      IF (LRND) 630,640,660                                             DVDQ0849
C     EXTREME ROUND-OFF ERROR--REDUCE E2H                               DVDQ0850
  630 K=(KBIT2/KQQ)-4                                                   DVDQ0851
      IF (K.LE.3) GO TO 640                                             DVDQ0852
      IF (K.GT.KQMAX) K=KQMAX                                           DVDQ0853
      E2H=E2H/PT(K+1)                                                   DVDQ0854
      GO TO 650                                                         DVDQ0855
  640 E2H=AMIN1(E2H,E2H*P3E1*E2HFAC)                                    DVDQ0856
  650 E2H=E2H*P1                                                        DVDQ0857
      TPS6=PTS4                                                         DVDQ0858
      GO TO 670                                                         DVDQ0859
C                                                                       DVDQ0860
  660 E2H=E2H*E2HFAC                                                    DVDQ0861
      TPS6=FLOAT(KQQ+2)                                                 DVDQ0862
C     TEST TO SEE IF DIFFERENCES DECREASE MORE RAPIDLY THAN NECESSARY   DVDQ0863
C                                                                       DVDQ0864
  670 IF (TPS5.LT.(TPS3*TPS6)) GO TO 680                                DVDQ0865
      IF (TPS2.LE.(TPS4*TPS6)) GO TO 760                                DVDQ0866
C     THEY DO   INCREASE KQ(I)                                          DVDQ0867
      IF (KQQ.NE.KQMAX) KQ(I)=KQ1                                       DVDQ0868
      GO TO 760                                                         DVDQ0869
C                                                                       DVDQ0870
C     TEST TO SEE IF DIFFERENCES DECREASE TOO SLOWLY                    DVDQ0871
  680 TPS6=TPS6*P25                                                     DVDQ0872
      IF ((TPS1.GT.(TPS3*TPS6)).OR.(TPS2.GT.(TPS4*TPS6))) GO TO 760     DVDQ0873
C     THEY DO                                                           DVDQ0874
      IF (LSTC.LE.0) GO TO 750                                          DVDQ0875
      IF (E2H.LT.P01) GO TO 750                                         DVDQ0876
      IF (LSC-LSTC) 690,750,770                                         DVDQ0877
  690 IF (KSTEP-4) 750,700,710                                          DVDQ0878
  700 KQ1=LSTC                                                          DVDQ0879
  710 LSC=KQ1                                                           DVDQ0880
C     END OF ONE DERIVATIVE EVALUATION PER STEP                         DVDQ0881
      GO TO 770                                                         DVDQ0882
C                                                                       DVDQ0883
C     AFTER HALVING H. REDUCE KQ(I) IF A DISCONTINUITY HAS OCCURRED.    DVDQ0884
  720 IF (LDOUB.EQ.(-2)) GO TO 340                                      DVDQ0885
      DT(KQQ+1,I)=D(KQQ+1)                                              DVDQ0886
      IF (LDOUB.EQ.(-1)) DT(KQQ+1,I)=D(KQQ+2)                           DVDQ0887
      K=KQQ                                                             DVDQ0888
  730 IF (K.EQ.1) GO TO 740                                             DVDQ0889
      IF ((DABS(D(K-1)).GT.(PT(2)*DABS(D(K+1)))).OR.                    DVDQ0890
     1  (DABS(D(K)).GT.(PT(2)*DABS(D(K+2))))) GO TO 740                 DVDQ0891
      K=K-1                                                             DVDQ0892
      GO TO 730                                                         DVDQ0893
  740 IF ((K+K).GE.KQQ) GO TO 780                                       DVDQ0894
      LDOUB=-4                                                          DVDQ0895
      E2H=0.                                                            DVDQ0896
      KQQ=K+1                                                           DVDQ0897
C                                                                       DVDQ0898
C                                                                       DVDQ0899
C     DIFFERENCES DECREASE TOO SLOWLY REDUCE KQ(I).                     DVDQ0900
  750 KQ(I)=KQQ-1                                                       DVDQ0901
      IF (KQQ.EQ.2) DT(3,I)=0.D0                                        DVDQ0902
  760 IF (E2H.LT.P01) GO TO 780                                         DVDQ0903
  770 LDOUB=2                                                           DVDQ0904
  780 CONTINUE                                                          DVDQ0905
C                                                                       DVDQ0906
C     IF THE OUTPUT OPTION IS ELIMINATED, REMOVE THE 6 FOLLOWING CARDS. DVDQ0907
      IF (NEQ.GT.0) GO TO 790                                           DVDQ0908
      IO2=MAX0(1,(KQQ-1))                                               DVDQ0909
      IO3=IO2+3                                                         DVDQ0910
      WRITE (6,5021) I,KQQ,LRND,LDOUB,E,E2H,EPS,                        DVDQ0911
     1  (IO1,D(IO1),IO1=IO2,IO3)                                        DVDQ0912
 5021 FORMAT (1H I2,I4,2I5,1PE13.3,2E11.3,4(3H  (,I2,1H),D10.3))        DVDQ0913
C                                                                       DVDQ0914
  790 CONTINUE                                                          DVDQ0915
C                                                                       DVDQ0916
C     END OF LOOP FOR CORRECTING, ESTIMATING THE ERROR, ETC.            DVDQ0917
C                                                                       DVDQ0918
C                                                                       DVDQ0919
C     IF THE INTERPOLATION CAPABILITY IS ELIMINATED REMOVE THE          DVDQ0920
C     FOLLOWING CARD.                                                   DVDQ0921
      IF (IFL.LT.0) GO TO 1250                                          DVDQ0922
C     TEST FOR HALVING H                                                DVDQ0923
      IF (LDOUB) 800,950,870                                            DVDQ0924
  800 LDOUB=LDOUB+1                                                     DVDQ0925
      IF (LDOUB+1) 810,870,820                                          DVDQ0926
  810 IF (LDOUB.EQ.(-2)) GO TO 820                                      DVDQ0927
C     ORDER IN AT LEAST ONE COMPONENT HAS BEEN GREATLY REDUCED          DVDQ0928
      LDOUB=0                                                           DVDQ0929
      GO TO 220                                                         DVDQ0930
  820 DO 860 I=1,NE                                                     DVDQ0931
      KQQ=KQ(I)                                                         DVDQ0932
      TP=DT(KQQ+1,I)                                                    DVDQ0933
      IF (KQQ.LE.3) GO TO 860                                           DVDQ0934
      IF (LDOUB.NE.0) GO TO 840                                         DVDQ0935
      DO 830 K=3,KQQ                                                    DVDQ0936
C     SECOND MODIFICATION OF DIFFERENCE TABLE AFTER HALVING H           DVDQ0937
  830 DT(K,I)=DT(K,I)+ETA(KQQ-1,K-2)*TP                                 DVDQ0938
      GO TO 860                                                         DVDQ0939
  840 DO 850 K=2,KQQ                                                    DVDQ0940
C     FIRST MODIFICATION OF DIFFERENCE TABLE AFTER HALVING H            DVDQ0941
  850 DT(K,I)=DT(K,I)+ETA(K-1,KQQ-1)*TP                                 DVDQ0942
  860 CONTINUE                                                          DVDQ0943
      IFL=0                                                             DVDQ0944
      GO TO 240                                                         DVDQ0945
C                                                                       DVDQ0946
  870 IFL=2                                                             DVDQ0947
      IF (LSTC.LE.0) GO TO 300                                          DVDQ0948
      IF (2-LSTC) 880,900,940                                           DVDQ0949
  880 LSTC=LSTC-1                                                       DVDQ0950
      IF (LSTC.EQ.3) GO TO 890                                          DVDQ0951
      IF (LSC) 920,960,920                                              DVDQ0952
  890 IFL=1                                                             DVDQ0953
      GO TO 300                                                         DVDQ0954
  900 IF (LSC-2) 910,930,920                                            DVDQ0955
  910 LSTC=0                                                            DVDQ0956
  920 LDOUB=2                                                           DVDQ0957
      GO TO 60                                                          DVDQ0958
  930 LSTC=1                                                            DVDQ0959
      LSC=0                                                             DVDQ0960
      GO TO 60                                                          DVDQ0961
  940 IF (LSC) 300,60,300                                               DVDQ0962
C                                                                       DVDQ0963
C     HALVE H                                                           DVDQ0964
  950 HH=FAC(2)*HH                                                      DVDQ0965
      IF (LSTC.LT.2) GO TO 990                                          DVDQ0966
      ERND=P25*ERND                                                     DVDQ0967
C     IN LOOP TO FIND A NEW INITIAL STEPSIZE                            DVDQ0968
      IF (ERND.GE.P1) GO TO 950                                         DVDQ0969
      EMAX = ERND                                                               
      LSTC=4                                                            DVDQ0970
  960 LSC=4                                                             DVDQ0971
      DO 970 I=1,NE                                                     DVDQ0972
  970 KQ(I)=1                                                           DVDQ0973
      IF (LSTC-3) 890,890,1170                                          DVDQ0974
C                                                                       DVDQ0975
C     ENTRY AFTER IFLAG=7                                               DVDQ0976
  980 IF (LDOUB.EQ.0) GO TO 990                                         DVDQ0977
      LSC=1                                                             DVDQ0978
      LSTC=1                                                            DVDQ0979
      GO TO 140                                                         DVDQ0980
C     TEST TO SEE IF H IS TOO SMALL FOR HALVING                         DVDQ0981
  990 IF (ABS(SNGL(HH)).GE.HMINA) GO TO 1040                            DVDQ0982
      IF (IFL.EQ.7) GO TO 1010                                          DVDQ0983
 1000 IFL=7                                                             DVDQ0984
      GO TO 1020                                                        DVDQ0985
C                                                                       DVDQ0986
 1010 HH=HH+HH                                                          DVDQ0987
      IFL=2                                                             DVDQ0988
 1020 H=HH                                                              DVDQ0989
      GO TO 310                                                         DVDQ0990
C                                                                       DVDQ0991
C                                                                       DVDQ0992
C     ERROR CRITERIA PERMIT DOUBLING                                    DVDQ0993
 1030 HH=HH+HH                                                          DVDQ0994
      IF (LSTC.EQ.1) GO TO 1050                                         DVDQ0995
      LSC=-3                                                            DVDQ0996
 1040 LSTC=-1                                                           DVDQ0997
C                                                                       DVDQ0998
C     CHANGE THE STEPSIZE                                               DVDQ0999
 1050 DO 1160 I=1,NE                                                    DVDQ1000
      KQQ=KQ(I)                                                         DVDQ1001
      IF (KQQ.NE.1) GO TO 1070                                          DVDQ1002
      DT(6,I)=0.D0                                                      DVDQ1003
      D(3)=DT(3,I)*PT(2)                                                DVDQ1004
      IF (D(3).GT.PT(3)) LSC=-6                                         DVDQ1005
      IF (LDOUB.NE.0) GO TO 1060                                        DVDQ1006
      KQM=8                                                             DVDQ1007
      IF (D(3).GE.PT(5)) DT(7,I)=DT(7,I)*PT(2)                          DVDQ1008
      D(3)=D(3)/PT(3)                                                   DVDQ1009
 1060 DT(3,I)=D(3)                                                      DVDQ1010
      GO TO 1160                                                        DVDQ1011
C                                                                       DVDQ1012
C     BEGINNING OF LOOP FOR CHANGING DIFFERENCE TABLE TO                DVDQ1013
C     CORRESPOND TO NEW VALUE OF H                                      DVDQ1014
 1070 DO 1080 K=1,KQQ                                                   DVDQ1015
      D(K)=DT(K,I)/PT(K)                                                DVDQ1016
 1080 IF (LDOUB.EQ.0) D(K)=D(K)/PT(K)                                   DVDQ1017
      KQQ2=KQQ-2                                                        DVDQ1018
      IF (KQQ2) 1160,1140,1090                                          DVDQ1019
 1090 DO 1130 J=1,KQQ2                                                  DVDQ1020
      IF (LDOUB.NE.0) GO TO 1110                                        DVDQ1021
C                                                                       DVDQ1022
C     HALVE                                                             DVDQ1023
      K=KQQ                                                             DVDQ1024
 1100 D(K-1)=D(K-1)+D(K)                                                DVDQ1025
      K=K-1                                                             DVDQ1026
      IF (K+J-KQQ) 1130,1130,1100                                       DVDQ1027
C                                                                       DVDQ1028
C     DOUBLE                                                            DVDQ1029
 1110 DO 1120 K=J,KQQ2                                                  DVDQ1030
 1120 D(K+1)=D(K+1)-D(K+2)                                              DVDQ1031
 1130 CONTINUE                                                          DVDQ1032
C                                                                       DVDQ1033
 1140 DO 1150 K=2,KQQ                                                   DVDQ1034
      IF (LDOUB.NE.0) D(K)=D(K)*PT(K)                                   DVDQ1035
      DT(K,I)=D(K)*PT(K)                                                DVDQ1036
 1150 CONTINUE                                                          DVDQ1037
C     DIFFERENCE TABLE NOW CORRESPONDS TO NEW VALUE OF H                DVDQ1038
C                                                                       DVDQ1039
 1160 CONTINUE                                                          DVDQ1040
 1170 H=HH                                                              DVDQ1041
      IF (LDOUB.NE.0) GO TO 50                                          DVDQ1042
      LFD=1                                                             DVDQ1043
      IF (LSTC.GE.0) GO TO 220                                          DVDQ1044
      LDOUB=-3                                                          DVDQ1045
      LSC=LSTC-KQM                                                      DVDQ1046
      GO TO 220                                                         DVDQ1047
C     END OF CHANGING STEPSIZE                                          DVDQ1048
C                                                                       DVDQ1049
C                                                                       DVDQ1050
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C  IN COLUMN ONE   DVDQ1051
C     OF THE 2 FOLLOWING CARDS                                          DVDQ1052
C1180 IF (7-IFL) 1181,980,220                                           DVDQ1053
C1181 IF (IFL-8) 60,1200,60                                             DVDQ1054
C     AND THEN REMOVE THE 2 FOLLOWING CARDS.                            DVDQ1055
 1180 K=IFL-5                                                           DVDQ1056
      GO TO (220,980,1200,1570,1570,1720,1720,60,1480,1450,1630,1570), KDVDQ1057
C                                                                       DVDQ1058
C     ILLEGAL VALUE OF PARAMETER   INTEGRATION CAN NOT PROCEED          DVDQ1059
 1190 IFL=8                                                             DVDQ1060
      GO TO 310                                                         DVDQ1061
 1200 WRITE (6,4000)                                                    DVDQ1062
 4000 FORMAT (26HOIFLAG=8 IN CALL TO DVDQ1.)                            DVDQ1063
      STOP                                                              DVDQ1064
C                                                                       DVDQ1065
C                                                                       DVDQ1066
 1210 IF (T-TFINAL) 200,1190,200                                        DVDQ1067
C                                                                       DVDQ1068
C     IF ONE DOES NOT WANT THE INTERPOLATION FEATURE, REMOVE ALL CARDS  DVDQ1069
C     BELOW THIS POINT (EXCEPT FOR THE END STATEMENT), AND ADD THE      DVDQ1070
C     FIVE FOLLOWING STATEMENTS.                                        DVDQ1071
C1220 IFL=4                                                             DVDQ1072
C     IF (TPD1.GT.TPD) GO TO 1280                                       DVDQ1073
C     GO TO 310                                                         DVDQ1074
C1280 IFL=3                                                             DVDQ1075
C     GO TO 310                                                         DVDQ1076
C                                                                       DVDQ1077
 1220 IFL=4                                                             DVDQ1078
      IF (KSTEP.NE.0) GO TO 1270                                        DVDQ1079
      TPD2=TPD                                                          DVDQ1080
C     ESTIMATE ERROR WHEN EXTRAPOLATION FROM INITIAL POINT IS REQUESTED DVDQ1081
 1230 HH=HH*TPD1*.75D0                                                  DVDQ1082
C                                                                       DVDQ1083
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE FOLLOWING CARD.    DVDQ1084
      IFLS=IFL                                                          DVDQ1085
      IFL=-1                                                            DVDQ1086
      GO TO 230                                                         DVDQ1087
C                                                                       DVDQ1088
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE 4 FOLLOWING CARDS. DVDQ1089
 1240 IF ((LGSD.EQ.0).OR.(IFLS.NE.4)) GO TO 20                          DVDQ1090
      LGSE=-1                                                           DVDQ1091
      TPD=FAC(1)                                                        DVDQ1092
      GO TO 1820                                                        DVDQ1093
 1250 HH=H                                                              DVDQ1094
      IF (EMAX.LT.P01) GO TO 1260                                       DVDQ1095
C     ERROR IS TOO LARGE, REDUCE H AND REPEAT THE FIRST STEP            DVDQ1096
      IF (TPD1.LT.0.D0) GO TO 1190                                      DVDQ1097
      LDOUB=1                                                                   
      ERND=FAC(1)/TPD1                                                          
      ERND=ERND*ERND*P25                                                        
C     SET FLAG SO INTERPOLATION IS DONE                                         
      IFLAG = 3                                                                 
      GO TO 950                                                         DVDQ1101
C                                                                       DVDQ1102
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C IN COLUMN ONE    DVDQ1103
C     OF THE FOLLOWING CARD                                             DVDQ1104
C1260 IFL=4                                                             DVDQ1105
C     AND THEN REMOVE THE 2 FOLLOWING CARDS.                            DVDQ1106
 1260 IFL=IFLS                                                          DVDQ1107
      IF (IFL.NE.4) GO TO 1790                                          DVDQ1108
      TPD=TPD2                                                          DVDQ1109
      IFLAG=3                                                           DVDQ1110
 1270 IF (TPD1.GT.TPD) GO TO 1280                                       DVDQ1111
      T=TFINAL                                                          DVDQ1112
      TPD=TPD1                                                          DVDQ1113
      GO TO 1290                                                        DVDQ1114
 1280 T=TOUT                                                            DVDQ1115
      IFL=3                                                             DVDQ1116
 1290 IF ((TPD.EQ.0.D0).AND.(IFLAG.LE.2)) GO TO 310                     DVDQ1117
C                                                                       DVDQ1118
C     INTERPOLATE FOR OUTPUT                                            DVDQ1119
 1300 TP=TPD                                                            DVDQ1120
      D(2)=TP                                                           DVDQ1121
      KQQ2=0                                                            DVDQ1122
      KDC=0                                                             DVDQ1123
      D(1)=PT(1)                                                        DVDQ1124
      DD(1)=PT(1)                                                       DVDQ1125
      DO 1310 K=2,KQM                                                   DVDQ1126
      DD(1)=DD(1)+PT(1)                                                 DVDQ1127
      TP=TP+PT(1)                                                       DVDQ1128
 1310 D(K+1)=(D(K)*TP)/DD(1)                                            DVDQ1129
      GO TO 1350                                                        DVDQ1130
C                                                                       DVDQ1131
C     COMPUTE THE INTERPOLATING INTEGRATION COEFFICIENTS                DVDQ1132
 1320 KQQ2=1                                                            DVDQ1133
      L=KQM-KDC                                                         DVDQ1134
      KDC=KDC+1                                                         DVDQ1135
 1330 IF (L.LE.0) GO TO 1350                                            DVDQ1136
      TP=0.D0                                                           DVDQ1137
      K=L                                                               DVDQ1138
      J=L+KDC                                                           DVDQ1139
 1340 JS=J-K                                                            DVDQ1140
      TP=TP+GAS(K)*D(JS+1)                                              DVDQ1141
      K=K-1                                                             DVDQ1142
      IF (K.GT.0) GO TO 1340                                            DVDQ1143
      D(J)=TP                                                           DVDQ1144
C                                                                       DVDQ1145
C     D(J) IS THE INTEGRATION COEFFICIENT FOR THE INTERPOLATION WHICH   DVDQ1146
C     CORRESPONDS TO GAM(J-KDC,KDC).                                    DVDQ1147
C                                                                       DVDQ1148
      L=L-1                                                             DVDQ1149
      GO TO 1330                                                        DVDQ1150
C     END OF COMPUTING INTEGRATION COEFFICIENTS                         DVDQ1151
C                                                                       DVDQ1152
C     PERFORM THE PARTIAL STEP INTEGRATION                              DVDQ1153
 1350 J=0                                                               DVDQ1154
      DO 1420 I=1,NE                                                    DVDQ1155
      IF (KDS.LE.0) KDD=IABS(KD(I))                                     DVDQ1156
      IF (KDC.GT.KDD) GO TO 1410                                        DVDQ1157
      TP=0.D0                                                           DVDQ1158
      KQQ=KQ(I)+KQQ2                                                    DVDQ1159
 1360 L=KQQ-KDC                                                         DVDQ1160
      IF (L.LE.0) GO TO 1370                                            DVDQ1161
      TP=TP+D(KQQ)*DT(L,I)                                              DVDQ1162
      KQQ=KQQ-1                                                         DVDQ1163
      IF (KQQ) 1390,1390,1360                                           DVDQ1164
 1370 K=J+KDD                                                           DVDQ1165
      L=KDC                                                             DVDQ1166
 1380 L=L-1                                                             DVDQ1167
      IF (L.EQ.0) GO TO 1400                                            DVDQ1168
      TP=TP*HH+YN(K)*FAC(L)*TPD**L                                      DVDQ1169
      K=K-1                                                             DVDQ1170
      GO TO 1380                                                        DVDQ1171
 1390 F(I)=TP                                                           DVDQ1172
      GO TO 1420                                                        DVDQ1173
 1400 Y(K)=YN(K)+HH*TP                                                  DVDQ1174
 1410 J=J+KDD                                                           DVDQ1175
 1420 CONTINUE                                                          DVDQ1176
      IF (KDC.NE.KDMAX) GO TO 1320                                      DVDQ1177
C     END OF PARTIAL STEP INTEGRATION                                   DVDQ1178
C                                                                       DVDQ1179
C     IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C IN COLUMN ONE    DVDQ1180
C     OF THE FOLLOWING CARD                                             DVDQ1181
C     GO TO 310                                                         DVDQ1182
C     ALL STATEMENTS BELOW THIS POINT SHOULD THEN BE REMOVED (EXCEPT    DVDQ1183
C     FOR THE END STATEMENT)                                            DVDQ1184
      IF (LGSE) 1800,310,1810                                           DVDQ1185
C                                                                       DVDQ1186
C                                                                       DVDQ1187
C     SECTION FOR COMPUTING GSTOPS                                      DVDQ1188
C                                                                       DVDQ1189
      ENTRY DVDQG(NG,NSTOP,G,GT)                                        DVDQ1190
C                                                                       DVDQ1191
C     VARIABLES IN THE CALLING SEQUENCE HAVE THE FOLLOWING TYPES.       DVDQ1192
C      INTEGER NG,NSTOP                                                  DVDQ1193
C      DOUBLE PRECISION G(1),GT(1)                                       DVDQ1194
C                                                                       DVDQ1195
C     A GSTOP IS DEFINED AS A RETURN WHICH IS MADE TO THE USER WHEN A   DVDQ1196
C     USER SPECIFIED FUNCTION G PASSES THROUGH ZERO. THE USER MAY       DVDQ1197
C     SPECIFY ANY NUMBER OF FUNCTIONS G OF TWO TYPES. ZEROS OF THE FIRSTDVDQ1198
C     TYPE ARE LOCATED WITHOUT REQUIRING A DERIVATIVE EVALUATION        DVDQ1199
C     BEYOND THE ZERO. THIS TYPE OF GSTOP REQUIRES THAT G BE EVALUATED  DVDQ1200
C     BEFORE EACH DERIVATIVE EVALUATION. ZEROS OF THE SECOND TYPE ARE   DVDQ1201
C     LOCATED USING INTERPOLATION, WHICH IS MORE ACCURATE THAN THE      DVDQ1202
C     EXTRAPOLATION USED IN THE PRECEDING CASE AND ONLY REQUIRES ONE    DVDQ1203
C     EVALUATION OF G PER STEP. THUS ONE SHOULD USE THE SECOND TYPE OF  DVDQ1204
C     GSTOP IF POSSIBLE. USERS NOT USING THE GSTOP FEATURE NEED READ    DVDQ1205
C     NO FURTHER.                                                       DVDQ1206
C                                                                       DVDQ1207
C     DVDQG IS USED AS A SET UP CALL TO INDICATE A CHANGE IN THE NUMBER DVDQ1208
C     OR TYPES OF GSTOPS. DVDQG SHOULD BE CALLED JUST BEFORE OR JUST    DVDQ1209
C     AFTER CALLING DVDQ IF                                             DVDQ1210
C     1. ONE WANTS TO TEST FOR GSTOPS BEGINNING WITH THE FIRST STEP.    DVDQ1211
C     2. A JOB IS BEING RUN AFTER ANOTHER JOB THAT USES THE GSTOP       DVDQ1212
C        FEATURE. DVDQG MUST BE CALLED EVEN IF ALL THE VARIABLES IN     DVDQ1213
C        THE NEW JOB ARE THE SAME.                                      DVDQ1214
C     IN ADDITION DVDQG MAY BE CALLED AT ANY TIME IN THE INTEGRATION    DVDQ1215
C     TO CHANGE THE NUMBER OR TYPE OF GSTOPS.                           DVDQ1216
C                                                                       DVDQ1217
C     THE USAGE OF THE VARIABLES IS GIVEN BELOW.                        DVDQ1218
C                                                                       DVDQ1219
C     NG=   THE NUMBER OF COMPONENTS IN G TO BE EXAMINED FOR A ZERO.    DVDQ1220
C           IF DVDQG IS CALLED AFTER THE FIRST STEP OF THE INTEGRATION, DVDQ1221
C           THEN G IS EVALUATED FOR THE FIRST TIME AT THE END OF THE    DVDQ1222
C           NEXT STEP AND THUS A GSTOP IS NOT DETECTED IF G CHANGES     DVDQ1223
C           SIGN ON THE CURRENT STEP. IF IT IS IMPORTANT THAT G BE      DVDQ1224
C           EVALUATED IMMEDIATELY SET NG EQUAL TO THE NEGATIVE OF THE   DVDQ1225
C           NUMBER OF COMPONENTS TO BE TESTED FOR A ZERO. SETTING NG    DVDQ1226
C           LESS THAN ZERO WHEN CALLING DVDQG BEFORE THE FIRST STEP IS  DVDQ1227
C           NOT NECESSARY AND IS LIABLE TO BE DISASTEROUS. IF DVDQG IS  DVDQ1228
C           CALLED DURING THE INTEGRATION THE FOLLOWING STATEMENT SHOULDDVDQ1229
C           BE A GO TO (THE COMPUTED GO TO FOLLOWING THE CALL TO DVDQ1).DVDQ1230
C                                                                       DVDQ1231
C     NSTOP=THE NUMBER OF COMPONENTS OF G THAT MUST BE EXAMINED FOR     DVDQ1232
C           A ZERO BEFORE COMPUTING THE DERIVATIVES (FIRST TYPE OF      DVDQ1233
C           GSTOP). IF NSTOP.LT.0 OR NSTOP.GT.ABS(NG) IFLAG IS SET      DVDQ1234
C           EQUAL 8 AND AN IMMEDIATE RETURN IS MADE. IF NSTOP.GT.0,     DVDQ1235
C           G(1),G(2),...,G(NSTOP) ARE EXAMINED FOR A ZERO BEFORE EACH  DVDQ1236
C           DERIVATIVE EVALUATION, THE REMAINING COMPONENTS (IF ANY)    DVDQ1237
C           ARE EXAMINED AT THE END OF EACH STEP. WHEN A GSTOP IS FOUND DVDQ1238
C           THE SUBROUTINE SETS NSTOP EQUAL TO THE INDEX OF THE         DVDQ1239
C           COMPONENT OF G RESPONSIBLE FOR THE STOP.                    DVDQ1240
C                                                                       DVDQ1241
C     G=    A VECTOR CONTAINING THE CURRENT VALUES OF THE FUNCTIONS     DVDQ1242
C           WHOSE ZEROS ARE TO BE DETERMINED.                           DVDQ1243
C                                                                       DVDQ1244
C     GT=   A VECTOR WITH THE SAME DIMENSION AS G USED BY THE           DVDQ1245
C           SUBROUTINE FOR TEMPORARY STORAGE.                           DVDQ1246
C                                                                       DVDQ1247
C     RETURNS FROM CALLING DVDQ1 WITH IFLAG.GT.8 SHOULD BE INTERPETED   DVDQ1248
C     AS FOLLOWS. (WE USE NSTOPI TO DENOTE THE INITIAL VALUE OF NSTOP.) DVDQ1249
C     IFLAG                                                             DVDQ1250
C     = 9 COMPUTE G(NSTOPI+1),...,G(ABS(NG)) (THE COMPONENTS OF G WITH  DVDQ1251
C         ZEROS TO BE LOCATED USING INTERPOLATION). THEN CALL DVDQ1.    DVDQ1252
C         NO RETURN IS MADE WITH IFLAG=9 IF NSTOPI=ABS(NG).             DVDQ1253
C     =10 COMPUTE G(1),G(2),...,G(NSTOPI) (THE COMPONENTS OF G WITH     DVDQ1254
C         ZEROS TO BE LOCATED USING EXTRAPOLATION). THEN CALL DVDQ1.    DVDQ1255
C         NO RETURN IS MADE WITH IFLAG=10 IF NSTOPI=0.                  DVDQ1256
C     =11 G(NSTOP) IS APPROXIMATELY ZERO. IF THERE ARE NO               DVDQ1257
C         DISCONTINUITIES SIMPLY CALL DVDQ1 TO CONTINUE THE INTEGRATION.DVDQ1258
C     =12 G(NSTOP) CHANGES SIGN, BUT THERE IS DIFFICULTY IN CONVERGING  DVDQ1259
C         TO A ZERO. THE USER MAY WISH TO MAKE A SPECIAL CHECK TO BE    DVDQ1260
C         CERTAIN THAT EVERYTHING IS ALL RIGHT. TO CONTINUE THE         DVDQ1261
C         INTEGRATION CALL DVDQ1.                                       DVDQ1262
C                                                                       DVDQ1263
C      DOUBLE PRECISION RG                                               DVDQ1264
C      DOUBLE PRECISION GI                                               DVDQ1265
C      DIMENSION GI(2),RG(3)                                             DVDQ1266
C                                                                       DVDQ1267
C     INITIALIZE FOR GSTOPS                                             DVDQ1268
      NGA=IABS(NG)                                                      DVDQ1269
      LGSS=-NGA                                                         DVDQ1270
      LGSD=0                                                            DVDQ1271
      LGSE=0                                                            DVDQ1272
      IFLG=-20                                                          DVDQ1273
      IF (NG.GE.0) RETURN                                               DVDQ1274
      IFLG=-IFL                                                         DVDQ1275
 1430 LGSD=NSTOP                                                        DVDQ1276
      IF (LGSD) 1190,1450,1440                                          DVDQ1277
 1440 IFL=15                                                            DVDQ1278
      GO TO 1470                                                        DVDQ1279
C     ENTRY WITH IFL=15                                                 DVDQ1280
 1450 LGSS=0                                                            DVDQ1281
      IF (LGSD-NGA) 1460,1480,1190                                      DVDQ1282
 1460 LGSS=LGSD+1                                                       DVDQ1283
      IFL=14                                                            DVDQ1284
 1470 IFLAG=IFL-5                                                       DVDQ1285
      GO TO 315                                                         DVDQ1286
C     ENTRY WITH IFL=14                                                 DVDQ1287
 1480 DO 1490 I=1,NGA                                                   DVDQ1288
 1490 GT(I)=G(I)                                                        DVDQ1289
      GO TO 1730                                                        DVDQ1290
C     END OF INITIALIZATION FOR GSTOPS                                  DVDQ1291
C                                                                       DVDQ1292
C     ENTRY TO EVALUATE G AT THE END OF THE STEP                        DVDQ1293
 1500 LGSE=1                                                            DVDQ1294
 1510 IGK=LGSS                                                          DVDQ1295
      IFLG=0                                                            DVDQ1296
      IFL=9                                                             DVDQ1297
      GO TO 310                                                         DVDQ1298
C     ENTRY TO EVALUATE G BEFORE EVALUATING THE DERIVATIVES             DVDQ1299
 1520 IFLG=IFL                                                          DVDQ1300
      IFL=10                                                            DVDQ1301
 1530 IFLAG=10                                                          DVDQ1302
      IGKM=LGSD                                                         DVDQ1303
 1540 IGK=1                                                             DVDQ1304
 1550 GO TO 315                                                         DVDQ1305
 1560 IGK=IGK+1                                                         DVDQ1306
      IF (IGK.GT.IGKM) GO TO 1650                                       DVDQ1307
C     ENTRY WITH IFL=9,10, AND 17                                       DVDQ1308
C     TEST FOR G CHANGING SIGN                                          DVDQ1309
 1570 IF (G(IGK)*GT(IGK)) 1600,1580,1590                                DVDQ1310
 1580 IF (GT(IGK).NE.0.D0) GO TO 1600                                   DVDQ1311
      IF (TL.EQ.TG) GO TO 1560                                          DVDQ1312
 1590 IF (LGSE.GT.0) GT(IGK)=G(IGK)                                     DVDQ1313
      GO TO 1560                                                        DVDQ1314
C     G CHANGES SIGN -- PREPARE FOR ITERATION TO FIND ZERO              DVDQ1315
 1600 NSTOP=IGK                                                         DVDQ1316
      NSTOPI=IGK                                                        DVDQ1317
      IFLGS=IFL                                                         DVDQ1318
C     COMPUTE INITIAL VALUE FOR RG (=RATIO OF PARTIAL STEPSIZE WHERE    DVDQ1319
C     G IS KNOWN/THE INTEGRATION STEPSIZE)                              DVDQ1320
      IF (IFLG.EQ.0) GO TO 1610                                         DVDQ1321
      RG(3)=FAC(1)                                                      DVDQ1322
      RG(2)=0.D0                                                        DVDQ1323
      IF ((IFLG.EQ.2).AND.(IGK.LT.LGSS)) RG(2)=FAC(1)                   DVDQ1324
      GO TO 1620                                                        DVDQ1325
 1610 RG(3)=0.D0                                                        DVDQ1326
      RG(2)=-FAC(1)                                                     DVDQ1327
 1620 IF (LGSE.LT.0) RG(3)=TPD                                          DVDQ1328
      LGSE=-3                                                           DVDQ1329
      GI(2)=GT(IGK)                                                     DVDQ1330
      EPSGS=RND                                                         DVDQ1331
      IFL=16                                                            DVDQ1332
      K=1                                                               DVDQ1333
      GO TO 1640                                                        DVDQ1334
C     END OF PREPARATION TO BEGIN THE ITERATION                         DVDQ1335
C                                                                       DVDQ1336
C     ENTRY WITH IFL=16                                                 DVDQ1337
C     ITERATE TO FIND GSTOP                                             DVDQ1338
 1630 K=1                                                               DVDQ1339
      IF ((GI(2)*G(IGK)).GT.0.D0) K=2                                   DVDQ1340
      IF (DABS(GI(K)).GT.DABS(G(IGK))) GO TO 1640                       DVDQ1341
C     CONVERGENCE PROBLEMS                                              DVDQ1342
      LGSE=LGSE-1                                                       DVDQ1343
      IF (LGSE.EQ.(-5)) EPSGS=PTS1                                      DVDQ1344
      EPSGS=EPSGS*PTS4                                                  DVDQ1345
 1640 GI(K)=G(IGK)                                                      DVDQ1346
      RG(K)=RG(3)                                                       DVDQ1347
C     SECANT ITERATION (GIVES NEW PARTIAL STEPSIZE/H)                   DVDQ1348
      TPD=RG(1)-(GI(1)*(RG(2)-RG(1)))/(GI(2)-GI(1))                     DVDQ1349
      T=TL+TPD*HH                                                       DVDQ1350
C     TEST FOR CONVERGENCE OF ITERATION                                 DVDQ1351
      IF (DABS(TPD-RG(3)).LE.EPSGS) GO TO 1560                          DVDQ1352
      RG(3)=TPD                                                         DVDQ1353
      GO TO 1300                                                        DVDQ1354
 1650 IF (10-IFL) 1660,1700,100                                         DVDQ1355
 1660 IF (IGKM.NE.NGA) GO TO 1710                                       DVDQ1356
      IF (LGSE.GT.(-3)) GO TO 1690                                      DVDQ1357
      IF (LSTC.NE.4) GO TO 1670                                         DVDQ1358
C     ESTIMATE ERROR -- GSTOP IS THE RESULT OF EXTRAPOLATING FROM       DVDQ1359
C     THE INITIAL POINT                                                 DVDQ1360
      TPD1=TPD                                                          DVDQ1361
      RG(3)=TPD                                                         DVDQ1362
      GO TO 1230                                                        DVDQ1363
 1670 IFL=11                                                            DVDQ1364
      IF (LGSE.LT.(-4)) IFL=12                                          DVDQ1365
 1680 IFLAG=IFL                                                         DVDQ1366
C     TEST TO SEE IF GSTOP IS PRECEDED BY ANOTHER STOP                  DVDQ1367
      IF (((T-TOUT)*HH.LE.0.D0).AND.((T-TFINAL)*HH.LE.0.D0)) GO TO 1300 DVDQ1368
C     IT IS                                                             DVDQ1369
      RG(3)=TPD                                                         DVDQ1370
      IFLS=IFL                                                          DVDQ1371
      GO TO 200                                                         DVDQ1372
 1690 LGSE=1                                                            DVDQ1373
      IFL=IFLG                                                          DVDQ1374
      IF (IFL.LT.0) GO TO 20                                            DVDQ1375
 1700 IGKM=NGA                                                          DVDQ1376
      IFL=IFLG                                                          DVDQ1377
      GO TO 310                                                         DVDQ1378
 1710 IFL=17                                                            DVDQ1379
      IFLAG=9                                                           DVDQ1380
      IGKM=NGA                                                          DVDQ1381
      GO TO 315                                                         DVDQ1382
C     ENTRY WITH IFL=11 AND 12                                          DVDQ1383
C     SET PARAMETERS TO INDICATE A GSTOP HAS BEEN FOUND                 DVDQ1384
 1720 GT(NSTOPI)=0.D0                                                   DVDQ1385
 1730 LGSE=1                                                            DVDQ1386
      IGKM=NGA                                                          DVDQ1387
      TG=TL                                                             DVDQ1388
      IF (IFLG) 1740,1760,1770                                          DVDQ1389
 1740 IF (IFL.LT.13) GO TO 1750                                         DVDQ1390
      IF (IFLG.EQ.(-20)) GO TO 100                                      DVDQ1391
      IFL=-IFLG                                                         DVDQ1392
      GO TO 310                                                         DVDQ1393
 1750 HH=H                                                              DVDQ1394
      GO TO 200                                                         DVDQ1395
 1760 TPD=0.D0                                                          DVDQ1396
      T=TL                                                              DVDQ1397
      LGSE=-2                                                           DVDQ1398
      GO TO 1300                                                        DVDQ1399
 1770 IF (IFLG-3) 220,200,200                                           DVDQ1400
 1780 IF (LGSE.EQ.(-1)) GO TO 1790                                      DVDQ1401
      LGSE=-1                                                           DVDQ1402
      GO TO 1220                                                        DVDQ1403
 1790 TPD=RG(3)                                                         DVDQ1404
      T=TL+TPD*HH                                                       DVDQ1405
      IF (LGSE.NE.(-1)) GO TO 1670                                      DVDQ1406
      IFL=IFLS                                                          DVDQ1407
      LGSE=-3                                                           DVDQ1408
      GO TO 1680                                                        DVDQ1409
 1800 IF (LGSE+2) 1550,1500,310                                         DVDQ1410
 1810 IF (TPD.LE.0.D0) GO TO 310                                        DVDQ1411
      LGSE=-2                                                           DVDQ1412
 1820 IFLG=IFL                                                          DVDQ1413
      IFL=17                                                            DVDQ1414
      IFLAG=9                                                           DVDQ1415
      IF (LGSD) 1540,1540,1530                                          DVDQ1416
C     END OF SECTION FOR COMPUTING GSTOPS                               DVDQ1417
C                                                                       DVDQ1418
C                                                                       DVDQ1419
C     IN SOME APPLICATIONS, FOR EXAMPLE MULTIPLE QUADRATURE, MORE THAN  DVDQ1420
C     ONE INTEGRATION SUBROUTINE IS REQUIRED. THIS IS NOT NECESSARY IF  DVDQ1421
C     ALL OF THE VARIABLES ASSOCIATED WITH ONE INTEGRATION ARE SAVED    DVDQ1422
C     OUTSIDE OF THE INTEGRATOR WHILE DOING OTHER INTEGRATIONS, AND     DVDQ1423
C     THEN RESTORING THEM WHEN NEEDED. THESE VARIABLES CAN BE RESTORED  DVDQ1424
C     BY CALLING AN ENTRY WHICH CONTAINS ALL OF THE VARIABLES IN THE    DVDQ1425
C     CALLING SEQUENCE AND ALL OF THE VARIABLES THAT MUST BE SAVED      DVDQ1426
C     WHENEVER CHANGING TO A DIFFERENT INTEGRATION. (THIS ENTRY MUST    DVDQ1427
C     BE ADDED BY THE USER AND SHOULD BE FOLLOWED BY A RETURN STATEMENT.DVDQ1428
C     AFTER CALLING THIS ENTRY EITHER DVDQ OR DVDQ1 SHOULD BE CALLED    DVDQ1429
C     DEPENDING ON WHETHER THE INTEGRATION IS BEING STARTED OR NOT.)    DVDQ1430
C     THE VARIABLES WHICH MUST BE SAVED ARE:                            DVDQ1431
C     NE,NV,KDS,KDMAX,KSOUT,LDOUB,LFD,LSC,LSTC,IFL,IFLS,KQM   (INTEGERS)DVDQ1432
C     ERND,QDEC,ERRMX,E2HAVE,E2HFAC,E2HMAX,RNDC   (REAL)                DVDQ1433
C     HH,TOUT,TL   (DOUBLE PRECISION)                                   DVDQ1434
C     IF THE GSTOP FEATURE IS USED, THE EVALUATION OF G REQUIRES AN     DVDQ1435
C     INTEGRATION, AND THIS INTEGRATION MAY RESULT IN ANOTHER GSTOP,    DVDQ1436
C     THEN SOME ADDITIONAL VARIABLES MUST BE SAVED.                     DVDQ1437
C     IN MANY APPLICATIONS NE(=NEQ),NV(=SUM OF ORDERS OF THE            DVDQ1438
C     DIFFERENTIAL EQUATIONS),KDS(=KD),AND KDMAX(=MAXIMUM ORDER OF ANY  DVDQ1439
C     DIFFERENTIAL EQUATION) WILL BE THE SAME FOR EVERY INTEGRATION,    DVDQ1440
C     AND HENCE NEED NOT BE SAVED.                                      DVDQ1441
C                                                                       DVDQ1442
C     INSTRUCTIONS FOR MAKING CERTAIN CHANGES IN THIS SUBROUTINE ARE    DVDQ1443
C     GIVEN THROUGHOUT THE LISTING. TO FIND THESE INSTRUCTIONS, SEE     DVDQ1444
C     BELOW.                                                            DVDQ1445
C                                                                       DVDQ1446
C     TO ELIMINATE THE GSTOP CAPABILITY, SEE JUST BELOW CARDS SEQUENCED DVDQ1447
C     406,552,594,640,1050,1083,1088,1102, AND 1179.                    DVDQ1448
C     THIS MAKES THE SUBROUTINE SHORTER AND REDUCES OVERHEAD A LITTLE.  DVDQ1449
C                                                                       DVDQ1450
C     TO REMOVE THE INTERPOLATION CAPABILITY, SEE JUST BELOW CARDS      DVDQ1451
C     SEQUENCED 919 AND 1069.                                           DVDQ1452
C     THE GSTOP FEATURE MUST ALSO BE ELIMINATED SINCE IT REQUIRES THE   DVDQ1453
C     INTERPOLATION CAPABILITY. IF OUTPUT POINTS ARE NOT HIT EXACTLY    DVDQ1454
C     (THEY ARE HIT EXACTLY IF HMAXA.LE.ABS(DELT), AND INITIAL H=       DVDQ1455
C     DELT*(2**(-K)), K=0,1,2...), THEN IFLAG=3 ON THE FIRST STEP THAT  DVDQ1456
C     (T-TOUT)*H.GT.0 (SEE THE USAGE OF DELT). IFLAG IS SET EQUAL TO 4  DVDQ1457
C     ON THE LAST STEP THAT (T-TFINAL)*H.LE.0.                          DVDQ1458
C                                                                       DVDQ1459
C     THE OUTPUT OPTION GIVES OUTPUT OF VARIABLES USED IN THE           DVDQ1460
C     INTEGRATION ON EVERY STEP THAT NEQ.LE.0. (WHICH OF COURSE MUST    DVDQ1461
C     BE SET AFTER THE INITIAL CALL TO THE INTEGRATOR)  TO ELIMINATE    DVDQ1462
C     THIS OPTION, SEE JUST BELOW CARDS SEQUENCED 683 AND 906.          DVDQ1463
C                                                                       DVDQ1464
C     THE CHECK OPTION WHEN ADDED TO THE OUTPUT OPTION OUTPUTS EVERY    DVDQ1465
C     VARIABLE IN THE CALLING SEQUENCE JUST AFTER ENTERING AND JUST     DVDQ1466
C     BEFORE LEAVING THE INTEGRATOR WHEN NEQ=0. THIS OUTPUT IS          DVDQ1467
C     SOMETIMES USEFUL IN DEBUGGING A PROGRAM. TO INCLUDE THIS OPTION   DVDQ1468
C     SEE JUST BELOW CARDS SEQUENCED 462 AND 651.                       DVDQ1469
C                                                                       DVDQ1470
      END                                                               DVDQ1471
CIBFTC HHFTDP  LIST  FORWARD TRIANGULARIZATION BY H.H. TRANS.  11 MAR 69HTD00100
